Citation
High-Resolution Spectral Analysis: The Missing Data Case

Material Information

Title:
High-Resolution Spectral Analysis: The Missing Data Case
Creator:
WANG, YANWEI ( Author, Primary )
Copyright Date:
2008

Subjects

Subjects / Keywords:
Apes ( jstor )
Estimate reliability ( jstor )
Estimated cost to complete ( jstor )
Estimators ( jstor )
Matrices ( jstor )
Missing data ( jstor )
Photographs ( jstor )
Preliminary estimates ( jstor )
Signals ( jstor )
Software applications ( jstor )

Record Information

Source Institution:
University of Florida
Holding Location:
University of Florida
Rights Management:
Copyright Yanwei Wang. Permission granted to the University of Florida to digitize, archive and distribute this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
Embargo Date:
12/31/2009
Resource Identifier:
656215755 ( OCLC )

Downloads

This item is only available as the following downloads:


Full Text

PAGE 1

HIGH-RESOLUTION SPECTRAL ANAL YSIS: THE MISSING D A T A CASE By Y ANWEI W ANG A DISSER T A TION PRESENTED TO THE GRADUA TE SCHOOL OF THE UNIVERSITY OF FLORID A IN P AR TIAL FULFILLMENT OF THE REQUIREMENTS F OR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORID A 2004

PAGE 2

Cop yrigh t 2004 b y Y an w ei W ang

PAGE 3

T o m y family

PAGE 4

A CKNO WLEDGMENTS I express m y sincere gratitude to m y advisor, Prof. Jian Li, who has guided me through the c hallenging journey of do ctoral studies. She has giv en me the freedom to explore new ideas, alw a ys encouraged me to ac hiev e more, and help ed me in all p ossible w a ys. I am esp ecially grateful to Prof. P etre Stoica for sharing his v ast and in v aluable exp erience with me. He has advised me with helpful suggestions, directed me to references, and taugh t me man y useful tec hniques. F urthermore, I am also indebted to Dr. Thomas Marzetta for sharing his insigh tful ideas with me. I w ould lik e to con v ey m y appreciation to Prof. Mark Sheplak, Prof. T oshik azu Nishida, and Prof. William Hager for b eing m y sup ervisory committee mem b ers. I am also grateful for their commen ts and suggestions on m y w ork. I am deeply thankful to m y family , who ha v e alw a ys supp orted me with their lo v e and understanding. In particular, I wish to express m y sp ecial thanks to m y wife, Guang Y ang, who has b een a source of strength and inspiration to me o v er these y ears. iv

PAGE 5

T ABLE OF CONTENTS page A CKNO WLEDGMENTS : : : : : : : : : : : : : : : : : : : : : : : : : : : : : iv LIST OF T ABLES : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : viii LIST OF FIGURES : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : ix ABSTRA CT : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : xii CHAPTER 1 INTR ODUCTION : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 1 1.1 Complete-Data Sp ectral Estimation . . . . . . . . . . . . . . . . . . . 1 1.2 Missing-Data Sp ectral Estimation . . . . . . . . . . . . . . . . . . . . 2 1.3 Main Con tributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.4 Thesis Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.4.1 Chapter 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.4.2 Chapter 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.4.3 Chapter 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.4.4 Chapter 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.4.5 Chapter 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.4.6 Chapter 7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2 APES F OR COMPLETE D A T A SPECTRAL ESTIMA TION : : : : : : : 7 2.1 In tro duction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2 Problem F orm ulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.3 F orw ard-Only APES Estimator . . . . . . . . . . . . . . . . . . . . . 8 2.4 Tw o-Step Filtering-Based In terpretation . . . . . . . . . . . . . . . . 10 2.5 Maxim um Lik eliho o d Fitting In terpretation . . . . . . . . . . . . . . 10 2.6 F orw ard-Bac kw ard Av eraging . . . . . . . . . . . . . . . . . . . . . . 13 2.7 F ast Implemen tation . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3 ONE-DIMENSIONAL MISSING-D A T A APES VIA EXPECT A TION MAX-IMIZATION::::::::::::::::::::::::::::::::::18 3.1 In tro duction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.2 EM for Missing-Data Sp ectral Estimation . . . . . . . . . . . . . . . 19 3.3 MAPES-EM1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.3.1 Exp ectation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.3.2 Maximization . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.4 MAPES-EM2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.4.1 Exp ectation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 v

PAGE 6

3.4.2 Maximization . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.5 Asp ects of In terest . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.5.1 Some Insigh ts in to the MAPES-EM Algorithms . . . . . . . . 28 3.5.2 MAPES-EM1 v ersus MAPES-EM2 . . . . . . . . . . . . . . . 28 3.5.3 Missing-Sample Estimation . . . . . . . . . . . . . . . . . . . 29 3.5.4 Initialization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.5.5 Stopping Criterion . . . . . . . . . . . . . . . . . . . . . . . . 30 3.6 MAPES Compared with GAPES . . . . . . . . . . . . . . . . . . . . 30 3.7 Numerical Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.8 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.9 App endix: GAPES . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.9.1 In tro duction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.9.2 GAPES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 4 TW O-DIMENSIONAL APES F OR COMPLETE D A T A SPECTRAL ESTIMATION::::::::::::::::::::::::::::::::::::48 4.1 In tro duction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 4.2 2-D APES Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 4.3 2-D ML-Based APES . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 5 TW O-DIMENSIONAL MAPES VIA EXPECT A TION MAXIMIZA TION AND CYCLICMAXIMIZATION::::::::::::::::::::::::::53 5.1 In tro duction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 5.2 2-D MAPES via Exp ectation Maximization . . . . . . . . . . . . . . 54 5.2.1 2-D MAPES-EM1 . . . . . . . . . . . . . . . . . . . . . . . . . 54 5.2.2 2-D MAPES-EM2 . . . . . . . . . . . . . . . . . . . . . . . . . 57 5.3 2-D MAPES via Cyclic Maximization . . . . . . . . . . . . . . . . . . 61 5.4 MAPES-EM v ersus MAPES-CM . . . . . . . . . . . . . . . . . . . . 63 5.4.1 P erformance Analysis . . . . . . . . . . . . . . . . . . . . . . . 63 5.4.2 Computational Analysis . . . . . . . . . . . . . . . . . . . . . 64 5.5 Numerical Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 5.5.1 Con v ergence Sp eed . . . . . . . . . . . . . . . . . . . . . . . . 66 5.5.2 P erformance Study . . . . . . . . . . . . . . . . . . . . . . . . 68 5.5.3 Syn thetic Ap erture Radar (SAR) Imaging Applications . . . . 69 5.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.7 App endix: 2-D GAPES . . . . . . . . . . . . . . . . . . . . . . . . . . 77 6 RANK-DEFICIENT R CF SPECTRAL ESTIMA TOR : : : : : : : : : : : 81 6.1 In tro duction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 6.2 Problem F orm ulation and Preliminaries . . . . . . . . . . . . . . . . . 83 6.3 Rank-Decien t Robust Cap on Filter-Bank Sp ectral Estimator . . . . 85 6.3.1 Rank-Decien t Sample Co v ariance Matrix . . . . . . . . . . . 86 6.3.2 Robust Cap on Filter-Bank (R CF) Approac h . . . . . . . . . . 87 6.4 Numerical Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 6.4.1 1-D Complex Sp ectral Estimation . . . . . . . . . . . . . . . . 95 6.4.2 2-D Complex Sp ectral Estimation . . . . . . . . . . . . . . . . 99 6.4.3 Com bining MAPES and R CF . . . . . . . . . . . . . . . . . . 101 6.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 6.6 App endixes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 6.6.1 2-D Complex Sp ectral Estimation . . . . . . . . . . . . . . . . 107 6.6.2 Rank-Decien t Cap on Beamformer via Co v ariance Fitting . . 109 vi

PAGE 7

6.6.3 The Conjugate Symmetry of the F orw ard-Bac kw ard FIR Filter 111 6.6.4 On the In terpla y Bet w een Sampling Errors and Steering Errors 112 6.6.5 F orm ulations of NCCF and HDI . . . . . . . . . . . . . . . . . 112 7 CONCLUSIONS AND FUTURE RESEAR CH : : : : : : : : : : : : : : : : 115 7.1 Summary of Results . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 7.2 T opics for F uture Researc h . . . . . . . . . . . . . . . . . . . . . . . . 116 REFERENCES : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 117 BIOGRAPHICAL SKETCH : : : : : : : : : : : : : : : : : : : : : : : : : : : : 123 vii

PAGE 8

LIST OF T ABLES T able page 3.1 RMSE's of the magnitude estimates obtained via the WFFT, GAPES, and MAPES-EM sp ectral estimators. . . . . . . . . . . . . . . . . . . 34 3.2 RMSE's of the phase (radian) estimates obtained via the WFFT, GAPES, and MAPES-EM sp ectral estimators. . . . . . . . . . . . . . 34 5.1 Computational complexities of the WFFT, GAPES, MAPES-EM, and MAPES-CM sp ectral estimators. . . . . . . . . . . . . . . . . . . . . 69 5.2 RMSE's of the magnitude estimates obtained via the WFFT, GAPES, MAPES-EM, and MAPES-CM sp ectral estimators. . . . . . . . . . . 70 5.3 RMSE's of the phase (radian) estimates obtained via the WFFT, GAPES, MAPES-EM, and MAPES-CM sp ectral estimators. . . . . . 70 6.1 RMSE's of the mo dulus and the phase (radian) estimates obtained b y the rank-decien t R CF and APES sp ectral estimators in the rst 1-D example. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 6.2 RMSE's of the frequency , the mo dulus, and the phase (radian) estimates obtained b y the rank-decien t R CF and APES sp ectral estimators in the second 1-D example. . . . . . . . . . . . . . . . . . . . . . 99 viii

PAGE 9

LIST OF FIGURES Figure page 3.1 Mo dulus of the missing-data sp ectral estimates ( N = 128, 2 n = 0 : 01, 51 (40%) missing samples): (a) T rue sp ectrum, (b) Complete-data APES, (c) WFFT, (d) GAPES with M = 64 and = 10 2 , (e) MAPES-EM1 with M = 64 and = 10 3 , and (f ) MAPES-EM2 with M = 64 and = 10 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.2 Mo dulus of the missing-data sp ectral estimates obtained via the MAPESEM algorithms at dieren t iterations ( N = 128, 2 n = 0 : 01, 51 (40%) missing samples): (a) MAPES-EM1, and (b) MAPES-EM2. . . . . . 37 3.3 In terp olation of the missing samples ( N = 128, 2 n = 0 : 01, 51 (40%) missing samples): (a) Real part of the data in terp olated via MAPESEM1, (b) Imaginary part of the data in terp olated via MAPES-EM1, (c) Real part of the data in terp olated via MAPES-EM2, and (d) Imaginary part of the data in terp olated via MAPES-EM2. . . . . . . . . 38 3.4 Estimation of the rst three missing samples vs. frequency ( N = 128, 2 n = 0 : 01, 51 (40%) samples are missing, v ertical dotted lines indicate the true frequency lo cations of the sp ectral comp onen ts with the closely spaced lines indicating the con tin uous sp ectral comp onen t): (a) Real part and (b) Imaginary part of the missing samples estimated via MAPES-EM2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.5 Mo dulus of the missing-data sp ectral estimates ( N = 128, 2 n = 0 : 01, 77 (60%) missing samples) obtained via: (a) WFFT, (b) GAPES with M = 64 and = 10 2 , (c) MAPES-EM1 with M = 64 and = 10 3 , and (d) MAPES-EM2 with M = 64 and = 10 3 . . . . . . . . . . . 40 3.6 Mo dulus of the missing-data sp ectral estimates ( N = 128, 2 n = 0 : 1, 51 (40%) missing samples) obtained via: (a) Complete-data WFFT, (b) Complete-data APES, (c) WFFT, (d) GAPES with M = 64 and = 10 2 , (e) MAPES-EM1 with M = 64 and = 10 3 , and (f ) MAPES-EM2 with M = 64 and = 10 3 . . . . . . . . . . . . . . . . 41 3.7 RMSE's of the estimates of the rst sp ectral line lo cated at f 1 = 0 : 05 Hz obtained via MAPES-EM1: (a) amplitude, and (b) phase (radian). 42 5.1 Mo dulus of the missing-data sp ectral estimates ( N 1 = 128, N 2 = 1, 51 (40%) missing samples): (a) Complete-data APES, (b) WFFT, (c) GAPES with M 1 = 64 and M 2 = 1, (d) MAPES-CM with M 1 = 64 and M 2 = 1, (e) MAPES-EM1 with M 1 = 64 and M 2 = 1, and (f ) MAPES-EM2 with M 1 = 64 and M 2 = 1. . . . . . . . . . . . . . . . 72 ix

PAGE 10

5.2 Mo dulus of the missing-data sp ectral estimates obtained at dieren t iterations ( N 1 = 128, N 2 = 1, 51 (40%) missing samples): (a) MAPESCM, (b) MAPES-EM1, (c) MAPES-EM2 , and (d) GAPES. . . . . . 73 5.3 Mo dulus of the 2-D sp ectra: (a) T rue sp ectrum, (b) 2-D completedata WFFT, (c) 2-D complete-data APES with M 1 = M 2 = 8, (d) 2-D data missing pattern, the blac k strip es indicate missing samples, (e) 2-D WFFT, (f ) 2-D GAPES with M 1 = M 2 = 8, (g) 2-D MAPESEM1 with M 1 = M 2 = 8, (h) 2-D MAPES-EM2 with M 1 = M 2 = 8, and (i) 2-D MAPES-CM with M 1 = M 2 = 8. . . . . . . . . . . . . . 74 5.4 Mo dulus of the SAR images of the Slicy ob ject obtained from a 32 32 data matrix with missing samples: (a) Photograph of the ob ject (tak en at 45 azim uth angle), (b) 2-D WFFT for a 288 288 (not 32 32) data matrix, (c) 2-D complete-data APES with M 1 = M 2 = 16, (d) 2-D data missing pattern, the blac k strip es indicate missing samples, (e) 2-D WFFT, (f ) 2-D GAPES with M 1 = M 2 = 16, and (g) 2-D MAPES-CM with M 1 = M 2 = 16. . . . . . . . . . . . . . . . . . . . 75 5.5 Mo dulus of the SAR images of the Bac kho e data obtained from a 32 32 data matrix with missing samples: (a) 3-D CAD mo del of the target and the 3-D K-space data dome, (b) mo dulus of the 2-D FFT image formed at 0 elev ation with 6 GHz bandwidth and a 110 azim uth cut (not 32 32), (c) 2-D complete-data WFFT, (d) 2-D complete-data APES with a 2-D lter of size 16 16, (e) 2-D data missing pattern, the blac k strip es indicate missing samples, (f ) 2-D WFFT, and (g) 2-D MAPES-CM with a 2-D lter of size 16 16. . 76 6.1 Mo dulus of 1-D sp ectral estimates ( N = 64): (a) T rue sp ectrum, (b) WFFT, (c) Cap on with M = 32, (d) APES with M = 32, (e) fullrank NCCF with M = 32 and = 0 : 3, (f ) rank-decien t NCCF with M = 56 and = 2, (g) full-rank R CF with M = 32 and = 0 : 3, and (h) rank-decien t R CF with M = 56 and = 0 : 3. . . . . . . . . . . . 102 6.2 Mo dulus of 1-D sp ectral estimates ( N = 64): (a) T rue sp ectrum, (b) APES with M = 32, (c) rank-decien t NCCF with M = 56 and = 2, and (d) rank-decien t R CF with M = 56 and = 0 : 3. . . . . . . . . 104 6.3 Mo dulus of the SAR images of the Slicy ob ject from a 24 24 data matrix: (a) Photograph of the ob ject (tak en at 45 azim uth angle), (b) 2-D WFFT with 288 288 (not 24 24) data matrix, (c) 2-D FFT, (d) 2-D WFFT, (e) 2-D CAPON with M = M = 12, (f ) 2-D APES with M = M = 12, (g) 2-D full-rank R CF with M = M = 12 and = 2, (h) 2-D rank-decien t R CF with M = M = 16 and = 2, (i) 2-D rank-decien t NCCF with M = M = 16 and = 0 : 2, and (j) 2-D HDI with M = M = 16 and % = 0 : 05. . . . . . . . . . . . . . . . . . 105 6.4 Mo dulus of the SAR images of the Slicy ob ject obtained from a 32 32 data matrix with missing samples: (a) 2-D MAPES-CM with M 1 = M 2 = 16, (b) 2-D MAPES-CM follo w ed b y 2-D rank-decien t R CF with a 20 20 lter-bank and a unit radius spherical constrain t. . . 114 x

PAGE 11

Abstract of Dissertation Presen ted to the Graduate Sc ho ol of the Univ ersit y of Florida in P artial F ul¯llmen t of the Requiremen ts for the Degree of Do ctor of Philosoph y HIGH-RESOLUTION SPECTRAL ANAL YSIS: THE MISSING D A T A CASE By Y an w ei W ang Decem b er 2004 Chair: Jian Li Ma jor Departmen t: Electrical and Computer Engineering Sp ectral estimation is imp ortan t in man y ¯elds including astronom y , meteorology , seismology , comm unications, economics, sp eec h analysis, medical imaging, radar, sonar, and underw ater acoustics. Most existing sp ectral estimation algorithms are devised for uniformly sampled complete-data sequences. Ho w ev er, the sp ectral estimation for data sequences with missing samples is also imp ortan t in a n um b er of applications ranging from astronomical time-series analysis to syn thetic ap erture radar imaging with angular div ersit y . The classical approac hes to sp ectral analysis include the discrete F ourier transform (DFT) or fast F ourier transform (FFT) based metho ds. These metho ds ha v e b een widely used for sp ectral estimation tasks due to their robustness and high computational e±ciency . Ho w ev er, they su®er from lo w resolution and p o or accuracy problems. Man y adv anced sp ectral estimation metho ds ha v e also b een prop osed, including parametric and nonparametric adaptiv e ¯ltering based approac hes. In general, the nonparametric approac hes are less sensitiv e to data mismo delling than their parametric coun terparts. Moreo v er, the adaptiv e ¯lter-bank based nonparametric xi

PAGE 12

sp ectral estimators can pro vide high resolution, lo w sidelob es, and accurate sp ectral estimates while retaining the robust nature of the nonparametric metho ds. This researc h is aimed to in v estigate adaptiv e ¯lter-bank based nonparametric complex sp ectral estimation of data sequences (one-dimensional case) or matrices (t w o-dimensional case) with missing samples o ccurring in arbitrary patterns. The main results of this w ork can b e summarized in three parts. First, w e consider the one-dimensional (1-D) case. W e dev elop t w o missingdata amplitude and phase estimation (MAPES) algorithms b y using a \maxim um lik eliho o d (ML) ¯tting" criterion. Then w e use the w ell-kno wn exp ectation maximization (EM) metho d to solv e the so-obtained estimation problem iterativ ely . Through n umerical sim ulations, w e demonstrate the excellen t p erformance of the MAPES algorithms for missing-data sp ectral estimation and missing data restoration. Next, w e presen t t w o-dimensional (2-D) extensions of the MAPES-EM algorithms in tro duced in the 1-D case. Then w e dev elop a new MAPES algorithm, referred to as MAPES-CM, b y solving an ML problem iterativ ely via cyclic maximization (CM). W e ha v e compared MAPES-EM with MAPES-CM and ha v e sho wn that MAPES-CM allo ws signi¯can t computational sa vings compared with MAPESEM, whic h is esp ecially desirable for 2-D applications. Numerical examples ha v e b een pro vided to demonstrate the p erformance of the 2-D MAPES algorithms. Finally , b y observing that all the missing-data algorithms dev elop ed previously estimate the missing data samples, w e prop ose to ac hiev e b etter sp ectral estimation p erformance, for example, higher resolution than APES, based on the (complete) data in terp olated via MAPES. It is kno wn that the Cap on sp ectral estimator has b etter resolution (narro w er mainlob es) compared with APES. In this part of the w ork, w e dev elop a rank-de¯cien t robust Cap on ¯lter-bank sp ectral estimator, whic h can ac hiev e m uc h higher resolution than the existing nonparametric sp ectral estimators. xii

PAGE 13

CHAPTER 1 INTR ODUCTION Sp ectral estimation is imp ortan t in man y ¯elds including astronom y , meteorology , seismology , comm unications, economics, sp eec h analysis, medical imaging, radar, sonar, and underw ater acoustics. Most existing sp ectral estimation algorithms are devised for uniformly sampled complete-data sequences [5, 81 , 9, 49 , 64 , 56 ]. Ho w ev er, the sp ectral estimation for data sequences with missing samples is also imp ortan t in a wide range of applications. F or example, sensor failure or outliers can lead to missingdata problems [23 ]. In astronomical, meteorological, or satellite-based applications, w eather or other conditions ma y disturb sample taking sc hemes (e.g., measuremen ts are only a v ailable during nigh ttime for astronomical applications), whic h will result in missing or gapp ed data [68 , 1, 44, 62 , 53 , 51, 41 , 63 , 59 , 8, 14, 76 , 27 ]. In syn thetic ap erture radar (SAR) imaging, missing sample problems arise when the syn thetic ap erture is gapp ed to reduce the radar resources needed for the high resolution imaging of a scene [30 , 32 , 29 ]. F or foliage p enetrating radar systems, certain radar op erating frequency bands are under strong electromagnetic or radio frequency (RF) in terference [61 , 60 ] so that the corresp onding samples m ust b e discarded resulting in missing data. Similar problems arise in data fusion via ultra wideband coheren t pro cessing [12 ]. 1.1 Complete-Data Sp ectral Estimation F or complete-data sp ectral estimation, extensiv e w ork has already b een carried out in the literature (see, e.g., the follo wing references: [71 , 25 , 42 ]). The con v en tional discrete F ourier transform (DFT) or fast F ourier transform (FFT) based metho ds ha v e b een widely used for sp ectral estimation tasks due to their robustness and 1

PAGE 14

2 high computational e±ciency . Ho w ev er, they su®er from lo w resolution and p o or accuracy problems. Man y adv anced sp ectral estimation metho ds ha v e also b een prop osed, including parametric [9 , 64 , 37 ] and nonparametric adaptiv e ¯ltering based approac hes [10 , 36 ]. In general, the nonparametric approac hes are less sensitiv e to data mismo delling than their parametric coun terparts. Moreo v er, the adaptiv e ¯lter-bank based nonparametric sp ectral estimators can pro vide high resolution, lo w sidelob es, and accurate sp ectral estimates while retaining the robust nature of the nonparametric metho ds [67 , 34 ]. They include the amplitude and phase estimation (APES) metho d [36 ] and the Cap on sp ectral estimator [10 ]. Ho w ev er, the complete-data sp ectral estimation metho ds do not w ork w ell in the missing-data case when the missing data samples are simply set to zero [33 ]. F or the DFT based sp ectral estimators, setting the missing samples to zero corresp onds to m ultiplying the original data with a windo wing function that assumes a v alue of one whenev er a sample is a v ailable, and zero otherwise. In the frequency domain, the resulting sp ectrum is the con v olution b et w een the F ourier transform of the complete data and that of the windo wing function. Since the F ourier transform of the windo wing function t ypically has an underestimated mainlob e and an extended pattern of undesirable sidelob es, the resulting sp ectrum will b e p o orly estimated and con tain sev ere artifacts. F or the parametric and adaptiv e ¯ltering based approac hes, similar p erformance degradations will also o ccur. 1.2 Missing-Data Sp ectral Estimation F or missing-data sp ectral estimation, v arious tec hniques ha v e b een dev elop ed previously . The CLEAN algorithm [21 , 65 , 53 ] is used to estimate the sp ectrum b y decon v olving the missing-data DFT sp ectrum (the so-called \dirt y map") in to the true signal sp ectrum (the so-called \true clean map") and the F ourier transform of the windo wing function (the so-called \dirt y b eam") via an iterativ e approac h. Although the CLEAN algorithm w orks for b oth missing and irregularly sampled data sequences, it

PAGE 15

3 cannot resolv e closely spaced sp ectral lines, and hence it ma y not b e a suitable to ol for high-resolution sp ectral estimation. The m ulti-tap er metho ds [7 , 17 ] compute sp ectral estimates b y assuming certain quadratic functions of the a v ailable data samples. The co e±cien ts in the corresp onding quadratic functions are optimized according to certain criteria, but it app ears that this approac h cannot o v ercome the resolution limit of DFT. T o ac hiev e high resolution, sev eral parametric algorithms, e.g., those based on an autoregressiv e (AR) or autoregressiv e mo ving-a v erage (ARMA) mo dels, w ere used to handle the missing-data problem [50 , 54, 55 , 57 , 6 , 16 , 83, 45, 20 , 22 , 18 ]. Although these parametric metho ds can pro vide impro v ed sp ectral estimates, they are sensitiv e to mo del errors. Recen tly , sev eral nonparametric adaptiv e ¯ltering based tec hniques ha v e b een dev elop ed for the missing-data problem. In [68 ] and [30 ], an extension of the APES metho d to the case of gapp ed data, referred to as GAPES, w as in tro duced. GAPES iterativ ely in terp olates the missing data and estimates the sp ectrum. In [29 ], the APES and Cap on sp ectral estimators w ere applied to p erio dically gapp ed data to obtain the PG-APES and PG-CAPON algorithms. Ho w ev er, these adaptiv e ¯ltering based metho ds can only deal with missing data o ccurring in gaps and do not w ork w ell for the more general problem of missing data samples o ccurring in arbitrary patterns. 1.3 Main Con tributions This researc h is aimed to in v estigate adaptiv e ¯lter-bank based nonparametric complex sp ectral estimation of data sequences (one-dimensional case) or matrices (t w o-dimensional case) with missing samples o ccurring in arbitrary patterns. The main results of this w ork can b e summarized in three parts. First, w e consider the one-dimensional (1-D) case. W e dev elop t w o missingdata amplitude and phase estimation (MAPES) algorithms b y using a \maxim um lik eliho o d (ML) ¯tting" criterion. Then w e use the w ell-kno wn exp ectation maximization (EM) metho d to solv e the so-obtained estimation problem iterativ ely . Through

PAGE 16

4 n umerical sim ulations, w e demonstrate the excellen t p erformance of the MAPES algorithms for missing-data sp ectral estimation and missing data restoration. Next, w e presen t t w o-dimensional (2-D) extensions of the MAPES-EM algorithms in tro duced in the 1-D case. Then w e dev elop a new MAPES algorithm, referred to as MAPES-CM, b y solving an ML problem iterativ ely via cyclic maximization (CM). W e ha v e compared MAPES-EM with MAPES-CM and ha v e sho wn that MAPES-CM allo ws signi¯can t computational sa vings compared with MAPESEM, whic h is esp ecially desirable for 2-D applications. Numerical examples ha v e b een pro vided to demonstrate the p erformance of the 2-D MAPES algorithms. Finally , b y observing that all the missing-data algorithms dev elop ed previously estimate the missing data samples, w e prop ose to ac hiev e b etter sp ectral estimation p erformance, for example, higher resolution than APES, based on the (complete) data in terp olated via MAPES. It is kno wn that the Cap on sp ectral estimator has b etter resolution (narro w er mainlob es) compared with APES. In this part of the w ork, w e dev elop a rank-de¯cien t robust Cap on ¯lter-bank sp ectral estimator, whic h can ac hiev e m uc h higher resolution than the existing nonparametric sp ectral estimators. 1.4 Thesis Outline The outlines of the remaining c hapters are giv en b elo w. 1.4.1 Chapter 2 In this c hapter, w e in tro duce the APES sp ectral estimator for the completedata case. W e deriv e the APES estimator from b oth pure narro wband-¯lter design and maxim um lik eliho o d ¯tting considerations. The APES estimator is needed for the missing-data algorithms dev elop ed in Chapter 3.

PAGE 17

5 1.4.2 Chapter 3 In this c hapter, w e dev elop t w o missing-data APES (MAPES) algorithms b y using a \ML ¯tting" criterion as discussed in Chapter 2. Then w e use the w ellkno wn exp ectation maximization (EM) metho d to solv e the so-obtained estimation problem iterativ ely . W e also demonstrate the adv an tage of MAPES-EM o v er GAPES b y comparing their design approac hes. 1.4.3 Chapter 4 This is an in tro ductory c hapter. W e pro vide the 2-D extension of the APES estimator. 1.4.4 Chapter 5 Tw o-dimensional extensions of the MAPES-EM algorithms are dev elop ed. Ho w ev er, due to the high computational complexit y in v olv ed, the direct application of MAPES-EM to large data sets, e.g., t w o-dimensional data, is computationally prohibitiv e. T o reduce the computational complexit y , w e dev elop another MAPES algorithm, referred to as MAPES-CM, b y solving a \ML ¯tting" problem iterativ ely via cyclic maximization (CM). MAPES-EM and MAPES-CM p ossess similar sp ectral estimation p erformance, but the computational complexit y of the latter is m uc h lo w er than that of the former. 1.4.5 Chapter 6 W e presen t a new nonparametric FIR ¯ltering based complex sp ectral estimator obtained b y using a rank-de¯cien t robust Cap on ¯lter-bank (R CF) approac h. W e compare the new sp ectral estimator with other FIR ¯ltering based approac hes including Cap on, APES, NCCF and HDI. It is found that b y allo wing the sample co v ariance matrix to b e rank-de¯cien t, rank-de¯cien t R CF can ac hiev e m uc h higher resolution than the existing sp ectral estimators.

PAGE 18

6 1.4.6 Chapter 7 The main results of the researc h are summarized in this c hapter. T opics for future researc h are also pro vided.

PAGE 19

CHAPTER 2 APES F OR COMPLETE D A T A SPECTRAL ESTIMA TION 2.1 In tro duction Filter-bank approac hes are commonly used for sp ectral analysis. As nonparametric sp ectral estimators, they attempt to compute the sp ectral con ten t of a signal without using an y a priori mo del information or making an y explicit mo del assumption ab out the signal. F or an y of these approac hes, the k ey elemen t is to design narro wband ¯lters cen tered at the frequencies of in terest. In fact, the w ell-kno wn p erio dogram can b e in terpreted as suc h a sp ectral estimator with a data-indep enden t ¯lter-bank. In general, data-dep enden t (or data-adaptiv e) ¯lters outp erform their data-indep enden t coun terparts and are hence preferred in man y applications. A w ell-kno wn adaptiv e ¯lter-bank metho d is the Cap on sp ectral estimator [10 ]. More recen tly , Li and Stoica [36 ] devised another adaptiv e ¯lter-bank metho d with enhanced p erformance, whic h is referred to as the Amplitude and Phase EStimation (APES). APES surpasses its riv als in sev eral asp ects [70 , 34 ] and ¯nds applications in v arious ¯elds [48 , 19, 82 , 58, 35 , 68]. In this c hapter, w e deriv e the APES estimator from b oth pure narro wband¯lter design [69 ] and maxim um lik eliho o d ¯tting [36 ] considerations. These materials pa v e the ground for the missing-data algorithms w e will presen t later. The remainder of this c hapter is organized as follo ws. The problem form ulation is giv en in Section 2.2 and the forw ard-only APES ¯lter is presen ted in Section 2.3. Section 2.4 pro vides a t w o-step ¯ltering in terpretation of the APES estimator. A maxim um lik eliho o d ¯tting in terpretation of the APES estimator is also in tro duced in Section 2.5. Section 2.6 sho ws ho w the forw ard-bac kw ard a v eraging can b e used to impro v e the p erformance 7

PAGE 20

8 of the estimator. A brief discussion ab out the fast implemen tation of APES app ears in Section 2.7. 2.2 Problem F orm ulation Consider the problem of estimating the amplitude sp ectrum of a complexv alued uniformly sampled discrete-time signal f y n g N ¡ 1 n =0 . F or a frequency ! of in terest, the signal y n is mo deled as y n = ® ( ! ) e j ! n + e n ( ! ) ; n = 0 ; :::; N ¡ 1 ; ! 2 [0 ; 2 ¼ ) ; (2.1) where ® ( ! ) denotes the complex amplitude of the sin usoidal comp onen t at frequency ! , and e n ( ! ) denotes the residual term (assumed zero-mean) whic h includes the unmo deled noise and in terference from frequencies other than ! . The problem of in terest is to estimate ® ( ! ) from f y n g N ¡ 1 n =0 for an y giv en frequency ! . 2.3 F orw ard-Only APES Estimator Let h ( ! ) denote the impulse resp onse of an M -tap FIR ¯lter-bank h ( ! ) = [ h 0 ( ! ) h 1 ( ! ) ¢ ¢ ¢ h M ¡ 1 ( ! )] T ; (2.2) where ( ¢ ) T denotes the transp ose. Then the ¯lter output can b e written as h H ( ! ) ¹ y l , where ¹ y l = [ y l y l +1 ¢ ¢ ¢ y l + M ¡ 1 ] T ; l = 0 ; :::; L ¡ 1 (2.3) are the M £ 1 o v erlapping forw ard data sub v ectors (snapshots) and L = N ¡ M + 1. Here ( ¢ ) H denotes the conjugate transp ose. F or eac h ! of in terest, w e consider the follo wing design ob jectiv e: min ® ( ! ) ; h ( ! ) L ¡ 1 X l =0 ¯ ¯ h H ( ! ) ¹ y l ¡ ® ( ! ) e j ! l ¯ ¯ 2 sub ject to h H ( ! ) a ( ! ) = 1 ; (2.4) where a ( ! ) is an M £ 1 v ector giv en b y a ( ! ) = [1 e j ! ¢ ¢ ¢ e j ! ( M ¡ 1) ] T : (2.5) In the ab o v e approac h, the ¯lter-bank h ( ! ) is designed suc h that

PAGE 21

9 a) the ¯ltered sequence is as close to a sin usoidal signal as p ossible in a least squares (LS) sense, b) the complex sp ectrum ® ( ! ) is not distorted b y the ¯ltering. Let ¹ g ( ! ) denote the normalized F ourier transform of ¹ y l : ¹ g ( ! ) = 1 L L ¡ 1 X l =0 ¹ y l e ¡ j ! l (2.6) and de¯ne ^ R = 1 L L ¡ 1 X l =0 ¹ y l ¹ y H l : (2.7) A straigh tforw ard calculation sho ws that the ob jectiv e function in (2.4) can b e rewritten as 1 L L ¡ 1 X l =0 ¯ ¯ h H ( ! ) ¹ y l ¡ ® ( ! ) e j ! l ¯ ¯ 2 = h H ( ! ) ^ R h ( ! ) ¡ ® ¤ ( ! ) h H ( ! ) ¹ g ( ! ) ¡ ® ( ! ) ¹ g H ( ! ) h ( ! ) + j ® ( ! ) j 2 = j ® ( ! ) ¡ h H ( ! ) ¹ g ( ! ) j 2 + h H ( ! ) ^ R h ( ! ) ¡ j h H ( ! ) ¹ g ( ! ) j 2 ; (2.8) where ( ¢ ) ¤ denotes the complex conjugate. The minimization of (2.8) with resp ect to ® ( ! ) is giv en b y ^ ® ( ! ) = h H ( ! ) ¹ g ( ! ) : (2.9) Insertion of (2.9) in (2.8) yields the follo wing minimization problem for the determination of h ( ! ): min h ( ! ) h H ( ! ) ^ S ( ! ) h ( ! ) sub ject to h H ( ! ) a ( ! ) = 1 (2.10) where ^ S ( ! ) 4 = ^ R ¡ ¹ g ( ! ) ¹ g H ( ! ) : (2.11) The solution to (2.10) is readily obtained [77 ] as ^ h ( ! ) = ^ S ¡ 1 ( ! ) a ( ! ) a H ( ! ) ^ S ¡ 1 ( ! ) a ( ! ) : (2.12)

PAGE 22

10 This is the forw ard-only APES ¯lter, and the forw ard-only APES estimator in (2.9) b ecomes ^ ® ( ! ) = a H ( ! ) ^ S ¡ 1 ( ! ) ¹ g ( ! ) a H ( ! ) ^ S ¡ 1 ( ! ) a ( ! ) : (2.13) 2.4 Tw o-Step Filtering-Based In terpretation The APES sp ectral estimator has a t w o-step ¯ltering in terpretation: passing the data f y n g N ¡ 1 n =0 through a bank of FIR bandpass ¯lters with v arying cen ter frequency ! , and then obtaining the sp ectrum estimate ^ ® ( ! ) for ! 2 [0 ; 2 ¼ ) from the ¯ltered data. F or eac h frequency ! , the corresp onding M -tap FIR ¯lter-bank is giv en b y (2.12). Hence the output obtained b y passing ¹ y l through the FIR ¯lter ^ h ( ! ) can b e written as ^ h H ( ! ) ¹ y l = ® ( ! )[ ^ h H ( ! ) a ( ! )] e j ! l + ¹ w l ( ! ) = ® ( ! ) e j ! l + ¹ w l ( ! ) ; (2.14) where ¹ w l ( ! ) = ^ h H ( ! ) ¹ e l ( ! ) denotes the residue term at the ¯lter output and the second equalit y follo ws from the iden tit y ^ h H ( ! ) a ( ! ) = 1 : (2.15) Th us, from the output of the FIR ¯lter, w e can obtain the least squares (LS) estimate of ® ( ! ) as ^ ® ( ! ) = ^ h H ( ! ) ¹ g ( ! ) : (2.16) 2.5 Maxim um Lik eliho o d Fitting In terpretation In this section, w e review the APES algorithm for complete-data sp ectral estimation follo wing the deriv ations in [36 ], whic h pro vide a \maxim um lik eliho o d (ML) ¯tting" in terpretation of the APES estimator.

PAGE 23

11 Recall the problem of estimating the amplitude sp ectrum of a complex-v alued uniformly sampled data sequence in tro duced in Section 2.2. The APES algorithm deriv ed b elo w estimates ® ( ! ) from f y n g N ¡ 1 n =0 for an y giv en frequency ! . P artition the data v ector y = [ y 0 y 1 ¢ ¢ ¢ y N ¡ 1 ] T (2.17) in to L o v erlapping sub v ectors (data snapshots) of size M £ 1 with the follo wing shifted structure ¹ y l = [ y l y l +1 ¢ ¢ ¢ y l + M ¡ 1 ] T ; l = 0 ; :::; L ¡ 1 ; (2.18) where L = N ¡ M + 1. Then, according to the data mo del in (2.1), the l th data snapshot ¹ y l can b e written as ¹ y l = ® ( ! ) a ( ! ) ¢ e j ! l + ¹ e l ( ! ) ; (2.19) where a ( ! ) is an M £ 1 v ector giv en b y (2.5) and ¹ e l ( ! ) = [ e l ( ! ) e l +1 ( ! ) ¢ ¢ ¢ e l + M ¡ 1 ( ! ) ] T . The APES algorithm mimics a ML approac h to estimate ® ( ! ) b y assuming that ¹ e l ( ! ), l = 0 ; 1 ; :::; L ¡ 1, are zero-mean circularly symmetric complex Gaussian random v ectors that are statistically indep enden t of eac h other and ha v e the same unkno wn co v ariance matrix Q ( ! ) = E £ ¹ e l ( ! ) ¹ e H l ( ! ) ¤ : (2.20) Then the co v ariance matrix of ¹ y l can b e written as R = j ® ( ! ) j 2 a ( ! ) a H ( ! ) + Q ( ! ) : (2.21) Since the v ectors f ¹ e l ( ! ) g L ¡ 1 l =0 in our case are o v erlapping, they are not statistically indep enden t of eac h other. Consequen tly , APES is not an exact ML estimator. Using the ab o v e assumptions, w e get the normalized surrogate log-lik eliho o d function of the data snapshots f ¹ y l g as follo ws: 1 L ln p ( f ¹ y l g j ® ( ! ) ; Q ( ! ) ) = ¡ M ln ¼ ¡ ln j Q ( ! ) j ¡ 1 L L ¡ 1 X l =0 [ ¹ y l ¡ ® ( ! )

PAGE 24

12 a ( ! ) e j ! l ¤ H Q ¡ 1 ( ! ) £ ¹ y l ¡ ® ( ! ) a ( ! ) e j ! l ¤ (2.22) = ¡ M ln ¼ ¡ ln j Q ( ! ) j ¡ tr ( Q ¡ 1 ( ! ) ¢ 1 L L ¡ 1 X l =0 £ ¹ y l ¡ ® ( ! ) a ( ! ) e j ! l ¤ £ ¹ y l ¡ ® ( ! ) a ( ! ) e j ! l ¤ H o ; (2.23) where tr f¢g and j ¢ j denote the trace and the determinan t of a matrix, resp ectiv ely . F or an y giv en ® ( ! ), maximizing (2.23) with resp ect to Q ( ! ) giv es ^ Q ® ( ! ) = 1 L L ¡ 1 X l =0 £ ¹ y l ¡ ® ( ! ) a ( ! ) e j ! l ¤ £ ¹ y l ¡ ® ( ! ) a ( ! ) e j ! l ¤ H : (2.24) Inserting (2.24) in to (2.23) yields the follo wing concen trated cost function (with c hanged sign) G = ¯ ¯ ¯ ^ Q ® ( ! ) ¯ ¯ ¯ = ¯ ¯ ¯ ¯ ¯ 1 L L ¡ 1 X l =0 £ ¹ y l ¡ ® ( ! ) a ( ! ) e j ! l ¤ £ ¹ y l ¡ ® ( ! ) a ( ! ) e j ! l ¤ H ¯ ¯ ¯ ¯ ¯ ; (2.25) whic h is to b e minimized with resp ect to ® ( ! ). By using the notation ¹ g ( ! ), ^ R , and ^ S ( ! ) de¯ned in (2.6), (2.7), and (2.11), resp ectiv ely , the cost function G in (2.25) b ecomes G = ¯ ¯ ¯ ^ R + j ® ( ! ) j 2 a ( ! ) a H ( ! ) ¡ ¹ g ( ! ) ® H ( ! ) a H ( ! ) ¡ ® ( ! ) a ( ! ) ¹ g H ( ! ) ¯ ¯ ¯ = ¯ ¯ ¯ ^ R ¡ ¹ g ( ! ) ¹ g H ( ! ) + [ ® ( ! ) a ( ! ) ¡ ¹ g ( ! )] [ ® ( ! ) a ( ! ) ¡ ¹ g ( ! )] H ¯ ¯ ¯ (2.26) = ¯ ¯ ¯ ^ S ( ! ) ¯ ¯ ¯ ¢ ¯ ¯ ¯ I + ^ S ¡ 1 ( ! ) [ ® ( ! ) a ( ! ) ¡ ¹ g ( ! )] [ ® ( ! ) a ( ! ) ¡ ¹ g ( ! )] H ¯ ¯ ¯ ; (2.27) where ^ S ( ! ) can b e recognized as an estimate of Q ( ! ). Making use of the iden tit y j I + AB j = j I + BA j , w e get G = ¯ ¯ ¯ ^ S ( ! ) ¯ ¯ ¯ ¢ n 1 + [ ® ( ! ) a ( ! ) ¡ ¹ g ( ! )] H ^ S ¡ 1 ( ! ) [ ® ( ! ) a ( ! ) ¡ ¹ g ( ! )] o : (2.28) Minimizing G with resp ect to ® ( ! ) yields ^ ® ( ! ) = a H ( ! ) ^ S ¡ 1 ( ! ) ¹ g ( ! ) a H ( ! ) ^ S ¡ 1 ( ! ) a ( ! ) : (2.29) Making use of the calculation in (2.26), w e get the estimate of Q ( ! ) as: ^ Q ( ! ) = ^ S ( ! ) + [ ^ ® ( ! ) a ( ! ) ¡ ¹ g ( ! )] [ ^ ® ( ! ) a ( ! ) ¡ ¹ g ( ! )] H : (2.30)

PAGE 25

13 In the APES algorithm, ^ ® ( ! ) is the sough t sp ectral estimate and ^ Q ( ! ) is the estimate of the n uisance matrix parameter Q ( ! ). The phrase \ML ¯tting criterion" used ab o v e can b e commen ted as follo ws. In some estimation problems, using the exact ML metho d is computationally prohibitiv e or ev en imp ossible. In suc h problems one can mak e a n um b er of simplifying assumptions and deriv e the corresp onding \ML" criterion. The estimates that minimize the so-obtained surrogate ML ¯tting criterion are not exact ML estimates, y et usually they ha v e go o d p erformance and generally they are b y design m uc h simpler to compute than the exact ML estimates. F or example, ev en if the data are not Gaussian distributed, a \ML" ¯tting criterion deriv ed under the Gaussian h yp othesis will often lead to computationally con v enien t and y et accurate estimates. Another example here is sin usoidal parameter estimation from data corrupted b y colored noise: the \ML" ¯tting criterion deriv ed under the assumption that the noise is white leads to parameter estimates of the sin usoidal comp onen ts whose accuracy asymptotically ac hiev es the exact CRB (deriv ed under the correct assumption of colored noise), see [73 , 66 ]. The APES metho d ([36 , 34]) is another example where a surrogate \ML" ¯tting criterion, deriv ed under the assumption that the data snapshots are Gaussian and indep enden t, leads to estimates with excellen t p erformance. W e follo w the same approac h in the follo wing c hapters b y extending the APES metho d to the missingdata case. 2.6 F orw ard-Bac kw ard Av eraging F orw ard-bac kw ard a v eraging has b een widely used for enhanced p erformance in man y sp ectral analysis applications. In the previous sections, w e obtained the APES ¯lter b y using only forw ard data v ectors. Here w e sho w that forw ard-bac kw ard a v eraging can b e readily incorp orated in to the APES ¯lter design b y considering b oth the forw ard and the bac kw ard data v ectors.

PAGE 26

14 Let the bac kw ard data sub v ectors (snapshots) b e constructed as ~ y l = [ y ¤ N ¡ l ¡ 1 y ¤ N ¡ l ¡ 2 ¢ ¢ ¢ y ¤ N ¡ l ¡ M ] T ; l = 0 ; :::; L ¡ 1 : (2.31) W e require that the outputs obtained b y running the data through the ¯lter b oth forw ard and bac kw ard are as close as p ossible to a sin usoid with frequency ! . This design ob jectiv e can b e written as min h ( ! ) ;® ( ! ) ;¯ ( ! ) 1 2 L L ¡ 1 X l =0 n ¯ ¯ h H ( ! ) ¹ y l ¡ ® ( ! ) e j ! l ¯ ¯ 2 + ¯ ¯ h H ( ! ) ~ y l ¡ ¯ ( ! ) e j ! l ¯ ¯ 2 o sub ject to h H ( ! ) a ( ! ) = 1 : (2.32) The minimization of (2.32) with resp ect to ® ( ! ) and ¯ ( ! ) giv es ^ ® ( ! ) = h H ( ! ) ¹ g ( ! ) and ^ ¯ ( ! ) = h H ( ! ) ~ g ( ! ), where ~ g ( ! ) is the normalized F ourier transform of ~ y l : ~ g ( ! ) = 1 L L ¡ 1 X l =0 ~ y l e ¡ j ! l : (2.33) It follo ws that (2.32) leads to min h ( ! ) h H ( ! ) ^ S f b ( ! ) h ( ! ) sub ject to h H ( ! ) a ( ! ) = 1 (2.34) where ^ S f b ( ! ) 4 = ^ R f b ¡ ¹ g ( ! ) ¹ g H ( ! ) + ~ g ( ! ) ~ g H ( ! ) 2 (2.35) with ^ R f = 1 L L ¡ 1 X l =0 ¹ y l ¹ y H l ; (2.36) ^ R b = 1 L L ¡ 1 X l =0 ~ y l ~ y H l ; (2.37) and ^ R f b = ^ R f + ^ R b 2 : (2.38) Note that here w e use ^ R f instead of ^ R to emphasize on the fact that it is estimated from the forw ard-only approac h. The solution of (2.34) is giv en b y ^ h f b ( ! ) = ^ S ¡ 1 f b ( ! ) a ( ! ) a H ( ! ) ^ S ¡ 1 f b ( ! ) a ( ! ) : (2.39)

PAGE 27

15 Due to the follo wing readily v eri¯ed relationship ~ y l = J ¹ y ¤ L ¡ l ¡ 1 ; (2.40) w e ha v e ~ g ( ! ) = J ¹ g ¤ ( ! ) ¢ e ¡ j ! ( L ¡ 1) ; (2.41) ^ R b = J ^ R T f J ; (2.42) and ~ g ( ! ) ~ g H ( ! ) = J [ ¹ g ( ! ) ¹ g H ( ! )] T J ; (2.43) where J denotes the exc hange matrix whose an tidiagonal elemen ts are ones and the remaining elemen ts are zeros. So ^ S f b ( ! ) can b e con v enien tly calculated as ^ S f b ( ! ) = ^ S f ( ! ) + J ^ S T f ( ! ) J 2 ; (2.44) where ^ S f ( ! ) 4 = ^ R f ¡ ¹ g ( ! ) ¹ g H ( ! ) : (2.45) Giv en the forw ard-bac kw ard APES ¯lter ^ h f b ( ! ), the forw ard-bac kw ard sp ectral estimator can b e written as ^ ® f b ( ! ) = a H ( ! ) ^ S ¡ 1 f b ( ! ) ¹ g ( ! ) a H ( ! ) ^ S ¡ 1 f b ( ! ) a ( ! ) : (2.46) Note that due to the ab o v e relationship, the forw ard-bac kw ard estimator of ¯ ( ! ) can b e simpli¯ed as ^ ¯ f b ( ! ) = h H f b ( ! ) ~ g ( ! ) = ^ ® ¤ f b ¢ e ¡ j ! ( N ¡ 1) ; (2.47) whic h indicates that from ^ ¯ f b ( ! ) w e will get the same forw ard-bac kw ard sp ectral estimator ^ ® f b ( ! ). In summary , the forw ard-bac kw ard APES ¯lter and APES sp ectral estimator still has the same forms as in (2.12) and (2.13), but ^ R and ^ S ( ! ) are replaced b y ^ R f b and ^ S f b ( ! ), resp ectiv ely . Note that ^ R f b and ^ S f b ( ! ) are p ersymmetric matrices.

PAGE 28

16 Compared with the non-p ersymmetric estimates ^ R f and ^ S f ( ! ), they are generally b etter estimates of the true R and Q ( ! ), where R and Q ( ! ) are the ideal co v ariance matrices with and without the presence of the signal of in terest, resp ectiv ely . F or simplicit y , all the APES-lik e algorithms w e dev elop in the subsequen t Chapters are based on the forw ard-only approac h. F or b etter estimation accuracy , the forw ard-bac kw ard a v eraging is used in all n umerical examples. 2.7 F ast Implemen tation The direct implemen tation of APES b y simply computing (2.13) for man y di®eren t ! of in terest is computationally demanding. Sev eral pap ers in the literature ha v e addressed this problem [39 , 82, 31 , 33 ]. Here w e giv e a brief discussion ab out implemen ting APES e±cien tly . T o a v oid the in v ersion of an M £ M matrix ^ S ( ! ) for eac h ! , w e use the matrix in v ersion lemma (see, e.g., [71 ]) to obtain ^ S ¡ 1 ( ! ) = ^ R ¡ 1 + ^ R ¡ 1 ¹ g ( ! ) ¹ g H ( ! ) ^ R ¡ 1 1 ¡ ¹ g H ( ! ) ^ R ¡ 1 ¹ g ( ! ) : (2.48) Let ^ R ¡ 1 = 2 denote the Cholesky factor of ^ R ¡ 1 , and let a R ( ! ) = ^ R ¡ 1 = 2 a ( ! ) g R ( ! ) = ^ R ¡ 1 = 2 ¹ g ( ! ) ³ ( ! ) = a H R ( ! ) a R ( ! ) % ( ! ) = a H R ( ! ) g R ( ! ) " ( ! ) = g H R ( ! ) g R ( ! ) : (2.49) Then w e can write (2.12) and (2.13) as ^ h ( ! ) = h ^ R ¡ 1 = 2 i H [(1 ¡ " ( ! )) a R ( ! ) + % ( ! ) g R ( ! )] ³ ( ! )(1 ¡ " ( ! )) + j % ( ! ) j 2 (2.50) and ^ ® ( ! ) = % ( ! ) ³ ( ! )(1 ¡ " ( ! )) + j % ( ! ) j 2 (2.51)

PAGE 29

17 whose implemen tation requires only the Cholesky factorization of the matrix ^ R that is indep enden t of ! . This strategy can b e readily generalized to the forw ard-bac kw ard a v eraging case. Since the complete-data case is not the fo cus of this b o ok, w e refer the readers to [39 , 82, 31 , 33] for more details ab out the e±cien t implemen tations of APES.

PAGE 30

CHAPTER 3 ONE-DIMENSIONAL MISSING-D A T A APES VIA EXPECT A TION MAXIMIZA TION 3.1 In tro duction The existing sp ectral estimation algorithms designed for uniformly sampled complete-data sequences p erform p o orly when applied to data sequences with missing samples if the missing samples are simply set to zero. Sev eral nonparametric algorithms ha v e recen tly b een dev elop ed to deal with the missing-data problem. They include, for example, GAPES [68 , 30] for gapp ed data and PG-APES [29 ] for p erio dically-gapp ed data. Ho w ev er, they are not really suitable for the general missing-data problem where the missing data samples o ccur in arbitrary patterns. F or example, GAPES can only deal with missing data o ccurring in gaps and PG-APES is designed for p erio dically-gapp ed data only . In this c hapter, w e consider the problem of nonparametric sp ectral estimation for data sequences with missing data samples o ccurring in arbitrary patterns (including the gapp ed-data case) [80 ]. W e dev elop t w o missing-data amplitude and phase estimation (MAPES) algorithms b y using a \maxim um lik eliho o d (ML) ¯tting" criterion as deriv ed in Chapter 2. Then w e use the w ell-kno wn exp ectation maximization (EM) [13 , 52 , 46 , 74 ] metho d to solv e the so-obtained estimation problem iterativ ely . Through n umerical sim ulations, w e demonstrate the excellen t p erformance of the MAPES algorithms for missing-data sp ectral estimation and missing data restoration. The remainder of this c hapter is organized as follo ws. In Section 3.2, w e giv e a brief review of the EM algorithm for the missing-data problem. In Section 3.3 and 18

PAGE 31

19 3.4, w e dev elop t w o nonparametric MAPES algorithms for the missing-data sp ectral estimation problem via the EM algorithm. Some asp ects of in terest are discussed in Section 3.5. In Section 3.6, w e compare MAPES with GAPES for the missing-data problem. Numerical results are pro vided in Section 3.7 to illustrate the p erformance of the MAPES-EM algorithms. Section 3.8 concludes this c hapter. 3.2 EM for Missing-Data Sp ectral Estimation Assume that some arbitrary samples of the uniformly sampled data sequence f y n g N ¡ 1 n =0 are missing. Due to these missing samples, whic h can b e treated as unkno wns, the surrogate log-lik eliho o d ¯tting criterion in (2.22) cannot b e maximized directly . W e sho w b elo w ho w to tac kle this general missing-data problem through the use of the EM algorithm. Recall that the g £ 1 v ector ° and the ( N ¡ g ) £ 1 v ector ¹ con tain all the a v ailable samples (incomplete data) and all the missing samples, resp ectiv ely , of the N £ 1 complete data v ector y . Then w e ha v e the follo wing relationships ° [ ¹ = f y n g N ¡ 1 n =0 (3.1) ° \ ¹ = Á; (3.2) where Á denotes the empt y set. Let µ = f ® ( ! ) ; Q ( ! ) g . An estimate ^ µ of µ can b e obtained b y maximizing the follo wing surrogate ML ¯tting criterion in v olving the a v ailable data v ector ° : ^ µ = arg max µ ln p ( ° j µ ) : (3.3) If ¹ w ere a v ailable, the ab o v e problem w ould b e easy to solv e (as sho wn in the previous section). In the absence of ¹ , ho w ev er, the EM algorithm maximizes the conditional (on ° ) exp ectation of the join t log-lik eliho o d function of ° and ¹ . The algorithm is iterativ e. A t the i th iteration, w e use ^ µ i ¡ 1 from the previous iteration to

PAGE 32

20 up date the parameter estimate b y maximizing the conditional exp ectation ^ µ i = arg max µ E n ln p ( ° ; ¹ j µ ) ¯ ¯ ¯ ° ; ^ µ i ¡ 1 o : (3.4) It can b e sho wn [74 , 11] that for eac h iteration, the increase in the surrogate loglik eliho o d function is greater than or equal to the increase in the exp ected join t surrogate log-lik eliho o d in (3.4), i.e., ln p ( ° j ^ µ i ) ¡ ln p ( ° j ^ µ i ¡ 1 ) ¸ E n ln p ( ° ; ¹ j ^ µ i ) ¯ ¯ ¯ ° ; ^ µ i ¡ 1 o ¡ E n ln p ( ° ; ¹ j ^ µ i ¡ 1 ) ¯ ¯ ¯ ° ; ^ µ i ¡ 1 o : (3.5) Since the data snapshots f ¹ y l g are o v erlapping, one missing sample ma y o ccur in man y snapshots (note that there is only one new sample b et w een t w o adjacen t data snapshots). So t w o approac hes are p ossible when w e try to estimate the missing data: estimate the missing data separately for eac h snapshot ¹ y l b y ignoring an y p ossible o v erlapping, or join tly for all snapshots f ¹ y l g L ¡ 1 l =0 b y observing the o v erlappings. In the follo wing t w o subsections, w e mak e use of these ideas to dev elop t w o di®eren t MAPES-EM algorithms, namely MAPES-EM1 and MAPES-EM2. 3.3 MAPES-EM1 In this section w e assume that the data snapshots f ¹ y l g L ¡ 1 l =1 are indep enden t of eac h other, and hence w e estimate the missing data separately for di®eren t data snapshots. F or eac h data snapshot ¹ y l , let ¹ ° l and ¹ ¹ l denote the v ectors con taining the a v ailable and missing elemen ts of ¹ y l , resp ectiv ely . In general, the indices of the missing comp onen ts could b e di®eren t for di®eren t l . Assume that ¹ ° l has dimension g l £ 1 where 1 · g l · M is the n um b er of a v ailable elemen ts in the snapshot ¹ y l . (Although g l could b e an y in teger that b elongs to the in terv al 0 · g l · M , w e assume for no w that g l 6 = 0. Later w e will explain what happ ens when g l = 0.) Then ¹ ° l and ¹ ¹ l are related to ¹ y l b y unitary transformations as follo ws: ¹ ° l = ¹ S T g ( l ) ¹ y l (3.6) ¹ ¹ l = ¹ S T m ( l ) ¹ y l ; (3.7)

PAGE 33

21 where ¹ S g ( l ) and ¹ S m ( l ) are M £ g l and M £ ( M ¡ g l ) unitary selection matrices suc h that ¹ S T g ( l ) ¹ S g ( l ) = I g l ; (3.8) ¹ S T m ( l ) ¹ S m ( l ) = I M ¡ g l ; (3.9) and ¹ S T g ( l ) ¹ S m ( l ) = 0 g l £ ( M ¡ g l ) : (3.10) F or example, if M = 5 and w e observ e the ¯rst, third, and fourth comp onen ts of ¹ y l , then g l = 3, ¹ S g ( l ) = 2 6 6 6 6 4 1 0 0 0 0 0 0 1 0 0 0 1 0 0 0 3 7 7 7 7 5 ; (3.11) and ¹ S m ( l ) = 2 6 6 6 6 4 0 0 1 0 0 0 0 0 0 1 3 7 7 7 7 5 : (3.12) Because w e clearly ha v e ¹ y l = £ ¹ S g ( l ) ¹ S T g ( l ) + ¹ S m ( l ) ¹ S T m ( l ) ¤ ¹ y l = ¹ S g ( l ) ¹ ° l + ¹ S m ( l ) ¹ ¹ l ; (3.13) the join t normalized surrogate log-lik eliho o d function of f ¹ ° l ; ¹ ¹ l g is obtained b y substituting (3.13) in to (2.23) 1 L ln p ( f ¹ ° l ; ¹ ¹ l g j ® ( ! ) ; Q ( ! ) ) = ¡ M ln ¼ ¡ ln j Q ( ! ) j ¡ tr ( Q ¡ 1 ( ! ) ¢ 1 L L ¡ 1 X l =0 £ ¹ S g ( l ) ¹ ° l + ¹ S m ( l ) ¹ ¹ l ¡ ® ( ! ) a ( ! ) e j ! l ¤ £ ¹ S g ( l ) ¹ ° l + ¹ S m ( l ) ¹ ¹ l ¡ ® ( ! ) a ( ! ) e j ! l ¤ H o : (3.14) Owing to the Gaussian assumption on ¹ y l , the follo wing random v ectors, · ¹ ¹ l ¹ ° l ¸ = · ¹ S T m ( l ) ¹ S T g ( l ) ¸ ¹ y l ; l = 0 ; :::; L ¡ 1 ; (3.15)

PAGE 34

22 are also Gaussian with mean · ¹ S T m ( l ) ¹ S T g ( l ) ¸ a ( ! ) ® ( ! ) e j ! l ; l = 0 ; :::; L ¡ 1 ; (3.16) and co v ariance matrix · ¹ S T m ( l ) ¹ S T g ( l ) ¸ Q ( ! ) £ ¹ S m ( l ) ¹ S g ( l ) ¤ ; l = 0 ; :::; L ¡ 1 : (3.17) F rom the Gaussian distribution of · ¹ ¹ l ¹ ° l ¸ , it follo ws that the probabilit y densit y function of ¹ ¹ l conditioned on ¹ ° l (for giv en µ = ^ µ i ¡ 1 ) is complex Gaussian with mean ¹ b l and co v ariance matrix ¹ K l [26 ]: ¹ ¹ l j ¹ ° l ; ^ µ i ¡ 1 » C N ( ¹ b l ; ¹ K l ) ; (3.18) where ¹ b l = E n ¹ ¹ l ¯ ¯ ¯ ¹ ° l ; ^ µ i ¡ 1 o = ¹ S T m ( l ) a ( ! ) ^ ® i ¡ 1 ( ! ) e j ! l + ¹ S T m ( l ) ^ Q i ¡ 1 ( ! ) ¹ S g ( l ) h ¹ S T g ( l ) ^ Q i ¡ 1 ( ! ) ¹ S g ( l ) i ¡ 1 ¡ ¹ ° l ¡ ¹ S T g ( l ) a ( ! ) ^ ® i ¡ 1 ( ! ) e j ! l ¢ (3.19) and ¹ K l = co v n ¹ ¹ l ¯ ¯ ¯ ¹ ° l ; ^ µ i ¡ 1 o = ¹ S T m ( l ) ^ Q i ¡ 1 ( ! ) ¹ S m ( l ) ¡ ¹ S T m ( l ) ^ Q i ¡ 1 ( ! ) ¹ S g ( l ) h ¹ S T g ( l ) ^ Q i ¡ 1 ( ! ) ¹ S g ( l ) i ¡ 1 ¹ S T g ( l ) ^ Q i ¡ 1 ( ! ) ¹ S m ( l ) : (3.20) 3.3.1 Exp ectation W e ev aluate the conditional exp ectation of the surrogate log-lik eliho o d in (3.14) using (3.18)-(3.20), whic h is most easily done b y adding and subtracting the conditional mean ¹ b l from ¹ ¹ l in (3.14) as follo ws: £ ¹ S g ( l ) ¹ ° l + ¹ S m ( l ) ¹ ¹ l ¡ ® ( ! ) a ( ! ) e j ! l ¤ = £ ¹ S m ( l )( ¹ ¹ l ¡ ¹ b l ) ¤ + £ ¹ S g ( l ) ¹ ° l + ¹ S m ( l ) ¹ b l ¡ ® ( ! ) a ( ! ) e j ! l ¤ : (3.21)

PAGE 35

23 The cross-terms that result from the expansion of the quadratic term in (3.14) v anish when w e tak e the conditional exp ectation. Therefore the exp ectation step yields E ½ 1 L ln p ( f ¹ ° l ; ¹ ¹ l g j ® ( ! ) ; Q ( ! ) ) ¯ ¯ ¯ f ¹ ° l g ; ^ ® i ¡ 1 ( ! ) ; ^ Q i ¡ 1 ( ! ) ¾ = ¡ M ln ¼ ¡ ln j Q ( ! ) j ¡ tr ( Q ¡ 1 ( ! ) 1 L L ¡ 1 X l =0 Ã ¹ S m ( l ) ¹ K l ¹ S T m ( l ) + £ ¹ S g ( l ) ¹ ° l + ¹ S m ( l ) ¹ b l ¡ ® ( ! ) a ( ! ) e j ! l ¤ £ ¹ S g ( l ) ¹ ° l + ¹ S m ( l ) ¹ b l ¡ ® ( ! ) a ( ! ) e j ! l ¤ H ´o : (3.22) 3.3.2 Maximization The maximization part of the EM algorithm pro duces up dated estimates for ® ( ! ) and Q ( ! ). The normalized exp ected surrogate log-lik eliho o d (3.22) can b e re-written as ¡ M ln ¼ ¡ ln j Q ( ! ) j ¡ tr n Q ¡ 1 ( ! ) 1 L P L ¡ 1 l =0 ³ ¹ ¡ l + £ ¹ z l ¡ ® ( ! ) a ( ! ) e j ! l ¤ £ ¹ z l ¡ ® ( ! ) a ( ! ) e j ! l ¤ H ´o ; (3.23) where w e ha v e de¯ned ¹ ¡ l 4 = ¹ S m ( l ) ¹ K l ¹ S T m ( l ) (3.24) and ¹ z l 4 = ¹ S g ( l ) ¹ ° l + ¹ S m ( l ) ¹ b l : (3.25) According to the deriv ation in Section 2.5, maximizing (3.23) with resp ect to ® ( ! ) and Q ( ! ) giv es ^ ® 1 ( ! ) = a H ( ! ) ¹ S ¡ 1 ( ! ) ¹ Z ( ! ) a H ( ! ) ¹ S ¡ 1 ( ! ) a ( ! ) (3.26) and ^ Q 1 ( ! ) = ¹ S ( ! ) + £ ^ ® 1 ( ! ) a ( ! ) ¡ ¹ Z ( ! ) ¤ £ ^ ® 1 ( ! ) a ( ! ) ¡ ¹ Z ( ! ) ¤ H ; (3.27) where ¹ Z ( ! ) 4 = 1 L L ¡ 1 X l =0 ¹ z l e ¡ j ! l (3.28) and ¹ S ( ! ) 4 = 1 L L ¡ 1 X l =0 ¹ ¡ l + 1 L L ¡ 1 X l =0 ¹ z l ¹ z H l ¡ ¹ Z ( ! ) ¹ Z H ( ! ) : (3.29)

PAGE 36

24 This completes the deriv ation of the MAPES-EM1 algorithm, a step-b y-step summary of whic h is giv en b elo w. Step 0 : Obtain an initial estimate of f ® ( ! ) ; Q ( ! ) g . Step 1 : Use the most recen t estimate of f ® ( ! ) ; Q ( ! ) g in (3.19) and (3.20) to calculate ¹ b l and ¹ K l , resp ectiv ely . Note that ¹ b l can b e regarded as the curren t estimate of the corresp onding missing samples. Step 2 : Up date the estimate of f ® ( ! ) ; Q ( ! ) g using (3.26) and (3.27). Step 3 : Rep eat Steps 1 and 2 un til practical con v ergence. Note that when g l = 0, whic h indicates that there is no a v ailable sample in the curren t data snapshot ¹ y l , ¹ S g ( l ) and ¹ ° l do not exist and ¹ S m ( l ) is an M £ M iden tit y matrix; hence, the ab o v e algorithm can still b e applied b y simply remo ving an y term that in v olv es ¹ S g ( l ) or ¹ ° l in the ab o v e equations. 3.4 MAPES-EM2 F ollo wing the observ ation that the same missing data ma y en ter in man y snapshots, w e prop ose a second metho d to implemen t the EM algorithm b y estimating the missing data sim ultaneously for all data snapshots. Recall that the a v ailable and missing data v ectors are denoted as ° ( g £ 1 v ector) and ¹ (( N ¡ g ) £ 1 v ector), resp ectiv ely . Let ¸ y denote the LM £ 1 v ector obtained b y concatenating all the snapshots ¸ y 4 = 2 6 4 ¹ y 0 . . . ¹ y L ¡ 1 3 7 5 = S g ° + S m ¹ ; (3.30) where S g ( LM £ g ) and S m ( LM £ ( N ¡ g )) are the corresp onding selection matrices for the a v ailable and missing data v ectors, resp ectiv ely . Due to the o v erlapping of f ¹ y l g , S g and S m are not unitary , but they are still orthogonal to eac h other: S T g S m = 0 g £ ( N ¡ g ) : (3.31)

PAGE 37

25 Instead of (3.6) and (3.7), w e ha v e from (3.30): ° = ( S T g S g ) ¡ 1 S T g ¸ y = ¸ S T g ¸ y (3.32) and ¹ = ( S T m S m ) ¡ 1 S T m ¸ y = ¸ S T m ¸ y : (3.33) The matrices ¸ S g and ¸ S m in tro duced ab o v e are de¯ned as ¸ S g 4 = S g ( S T g S g ) ¡ 1 (3.34) and ¸ S m 4 = S m ( S T m S m ) ¡ 1 ; (3.35) and they are also orthogonal to eac h other: ¸ S T g ¸ S m = 0 g £ ( N ¡ g ) : (3.36) Note that S T g S g and S T m S m are diagonal matrices where eac h diagonal elemen t indicates ho w man y times the corresp onding sample app ears in ¸ y due to the o v erlapping of f ¹ y l g . Hence b oth S T g S g and S T m S m can b e easily in v erted. No w the normalized surrogate log-lik eliho o d function in (2.22) can b e written as 1 L ln p ( ¸ y j ® ( ! ) ; Q ( ! ) ) = ¡ M ln ¼ ¡ 1 L ln j D ( ! ) j ¡ 1 L [ ¸ y ¡ ® ( ! ) ½ ( ! )] H D ¡ 1 ( ! ) [ ¸ y ¡ ® ( ! ) ½ ( ! )] ; (3.37) where ½ ( ! ) and D ( ! ) are de¯ned as: ½ ( ! ) 4 = 2 6 4 e j ! 0 a ( ! ) . . . e j ! ( L ¡ 1) a ( ! ) 3 7 5 (3.38) and D ( ! ) 4 = 2 6 4 Q ( ! ) 0 . . . 0 Q ( ! ) 3 7 5 : (3.39)

PAGE 38

26 Substituting (3.30) in to (3.37), w e obtain the join t surrogate log-lik eliho o d of ° and ¹ : 1 L ln p ( ° ; ¹ j ® ( ! ) ; Q ( ! ) ) = 1 L ½ ¡ LM ln ¼ ¡ ln j D ( ! ) j ¡ [ S g ° + S m ¹ ¡ ® ( ! ) ½ ( ! )] H D ¡ 1 ( ! ) [ S g ° + S m ¹ ¡ ® ( ! ) ½ ( ! )] o + C J ; (3.40) where C J is a constan t that accoun ts for the Jacobian of the non unitary transformation b et w een ¸ y and ° and ¹ in (3.30). T o deriv e the EM algorithm for the curren t set of assumptions, w e note that for giv en ^ ® i ¡ 1 ( ! ) and ^ Q i ¡ 1 ( ! ), w e ha v e (as in (3.18)-(3.20)): ¹ j ° ; ^ µ i ¡ 1 » C N ( b ; K ) (3.41) where b = E n ¹ ¯ ¯ ¯ ° ; ^ µ i ¡ 1 o = ¸ S T m ½ ( ! ) ^ ® i ¡ 1 ( ! ) + ¸ S T m ^ D i ¡ 1 ( ! ) ¸ S g h ¸ S T g ^ D i ¡ 1 ( ! ) ¸ S g i ¡ 1 ³ ° ¡ ¸ S T g ½ ( ! ) ^ ® i ¡ 1 ( ! ) ´ (3.42) and K = co v n ¹ ¯ ¯ ¯ ° ; ^ µ i ¡ 1 o = ¸ S T m ^ D i ¡ 1 ( ! ) ¸ S m ¡ ¸ S T m ^ D i ¡ 1 ( ! ) ¸ S g h ¸ S T g ^ D i ¡ 1 ( ! ) ¸ S g i ¡ 1 ¸ S T g ^ D i ¡ 1 ( ! ) ¸ S m : (3.43) 3.4.1 Exp ectation F ollo wing the same steps as in (3.21)-(3.22), w e obtain the conditional exp ectation of the surrogate log-lik eliho o d function in (3.40): E ½ 1 L ln p ( ° ; ¹ j ® ( ! ) ; Q ( ! ) ) ¯ ¯ ¯ ° ; ^ ® i ¡ 1 ( ! ) ; ^ Q i ¡ 1 ( ! ) ¾ = ¡ M ln ¼ ¡ 1 L ln j D ( ! ) j ¡ tr ½ 1 L D ¡ 1 ( ! ) µ S m KS T m + [ S g ° + S m b ¡ ® ( ! ) ½ ( ! )] [ S g ° + S m b ¡ ® ( ! ) ½ ( ! )] H ´o + C J : (3.44)

PAGE 39

27 3.4.2 Maximization T o maximize the exp ected surrogate log-lik eliho o d function in (3.44), w e need to exploit the kno wn structure of D ( ! ) and ½ ( ! ). Let 2 6 4 z 0 . . . z L ¡ 1 3 7 5 4 = S g ° + S m b (3.45) denote the data snapshots made up of the a v ailable and estimated data samples, where eac h z l , l = 0 ; :::; L ¡ 1, is an M £ 1 v ector. Also let ¡ 0 , ..., ¡ L ¡ 1 b e the M £ M blo c ks on the blo c k diagonal of S m KS T m . Then the exp ected surrogate loglik eliho o d function w e need to maximize with resp ect to ® ( ! ) and Q ( ! ) b ecomes (to within an additiv e constan t) ¡ ln j Q ( ! ) j ¡ tr ( Q ¡ 1 ( ! ) 1 L L ¡ 1 X l =0 ³ ¡ l + £ z l ¡ ® ( ! ) a ( ! ) e j ! l ¤ £ z l ¡ ® ( ! ) a ( ! ) e j ! l ¤ H ´ ) ; (3.46) The solution can b e readily obtained b y a deriv ation similar to that in Section 2.5: ^ ® 2 ( ! ) = a H ( ! ) ¸ S ¡ 1 ( ! ) Z ( ! ) a H ( ! ) ¸ S ¡ 1 ( ! ) a ( ! ) (3.47) and ^ Q 2 ( ! ) = ¸ S ( ! ) + [ ^ ® 2 ( ! ) a ( ! ) ¡ Z ( ! )] [ ^ ® 2 ( ! ) a ( ! ) ¡ Z ( ! )] H ; (3.48) where S ( ! ) and Z ( ! ) are de¯ned as ¸ S ( ! ) 4 = 1 L L ¡ 1 X l =0 ¡ l + 1 L L ¡ 1 X l =0 z l z H l ¡ Z ( ! ) Z H ( ! ) (3.49) and Z ( ! ) 4 = 1 L L ¡ 1 X l =0 z l e ¡ j ! l : (3.50) The deriv ation of the MAPES-EM2 algorithm is th us complete, and a step-b y-step summary of this algorithm follo ws. Step 0 : Obtain an initial estimate of f ® ( ! ) ; Q ( ! ) g .

PAGE 40

28 Step 1 : Use the most recen t estimates of f ® ( ! ) ; Q ( ! ) g in (3.42) and (3.43) to calculate b and K . Note that b can b e regarded as the curren t estimate of the missing sample v ector. Step 2 : Up date the estimates of f ® ( ! ) ; Q ( ! ) g using (3.47) and (3.48). Step 3 : Rep eat Steps 1 and 2 un til practical con v ergence. 3.5 Asp ects of In terest 3.5.1 Some Insigh ts in to the MAPES-EM Algorithms Comparing f ^ ® 1 ( ! ) ; ^ Q 1 ( ! ) g (or f ^ ® 2 ( ! ) ; ^ Q 2 ( ! ) g ) with f ^ ® ( ! ) ; ^ Q ( ! ) g , w e can see that the EM algorithms are doing some in tuitiv ely ob vious things. In particular, the estimator of ® ( ! ) estimates the missing data and then uses the estimate f ¹ b l g (or b ) as though it w ere correct. The estimator of Q ( ! ) do es the same thing, but it also adds an extra term in v olving the conditional co v ariance ¹ K l (or K ), whic h can b e regarded as a generalized diagonal loading op eration to mak e the sp ectral estimate robust against estimation errors. W e stress again that the MAPES approac h is based on a surrogate lik eliho o d function whic h is not the true lik eliho o d of the data snapshots. Ho w ev er, suc h surrogate lik eliho o d functions (for instance, based on false uncorrelatedness or Gaussian assumptions) are kno wn to lead to satisfactory ¯tting criteria, under fairly reasonable conditions (see, e.g., [40 , 73 ]). F urthermore, it can b e sho wn that the EM algorithm applied to suc h a surrogate lik eliho o d function (whic h is a v alid probabilit y distribution function) still has the k ey prop ert y in (3.5) to monotonically increase the function at eac h iteration. 3.5.2 MAPES-EM1 v ersus MAPES-EM2 Due to the fact that at eac h iteration and at eac h frequency of in terest ! , MAPES-EM2 estimates the missing samples only once (for all data snapshots), it has

PAGE 41

29 a lo w er computational complexit y than MAPES-EM1, whic h estimates the missing samples separately for eac h data snapshot. It is also in teresting to observ e that MAPES-EM1 mak es the assumption that the snapshots f ¹ y l g are indep enden t when form ulating the surrogate data lik eliho o d function, and it main tains this assumption when estimating the missing data | hence a \consisten t" ignoring of the o v erlapping. On the other hand, MAPES-EM2 mak es the same assumption when form ulating the surrogate data lik eliho o d function, but in a somewhat \inconsisten t" manner it observ es the o v erlapping when estimating the missing data. This suggests that MAPES-EM2, whic h estimates few er unkno wns than MAPES-EM1, ma y not necessarily ha v e a (m uc h) b etter p erformance, as migh t b e exp ected, see the examples in Section 3.7. 3.5.3 Missing-Sample Estimation F or man y applications, suc h as data restoration, estimating the missing samples is needed and can b e done via the MAPES-EM algorithms. F or MAPES-EM2, at eac h frequency of in terest ! , w e tak e the conditional mean b as an estimate of the missing sample v ector. The ¯nal estimate of the missing sample v ector is the a v erage of the b 's obtained from all frequencies of in terest. F or MAPES-EM1, at eac h frequency ! of in terest, there are m ultiple estimates (obtained from di®eren t o v erlapping data snapshots) for the same missing sample. W e calculate the mean of these m ultiple estimates b efore a v eraging once again across all frequencies of in terest. W e remark that w e should not consider the f ¹ b l g (or b ) at eac h frequency ! as an estimate of the ! -comp onen t of the missing data b ecause other frequency comp onen ts con tribute to the residue term as w ell, whic h determines the co v ariance matrix Q ( ! ) in the APES mo del.

PAGE 42

30 3.5.4 Initialization Since in general there is no guaran tee that the EM algorithm will con v erge to a global maxim um, the MAPES-EM algorithms ma y con v erge to a lo cal maxim um whic h dep ends on the initial estimate ^ µ 0 used. T o demonstrate the robustness of our MAPES-EM algorithms to the c hoice of the initial estimate, w e will simply let the initial estimate of ® ( ! ) b e giv en b y the windo w ed FFT (WFFT) with the missing data samples set to zero. The initial estimate of Q ( ! ) follo ws from (2.24), where again, the missing data samples are set to zero. 3.5.5 Stopping Criterion W e stop the iteration of the MAPES-EM algorithms whenev er the relativ e c hange in the total p o w er of the sp ectra corresp onding to the curren t ( ^ ® i ( ! k )) and previous ( ^ ® i ¡ 1 ( ! k )) estimates is smaller than a preselected threshold (e.g., ² = 10 ¡ 3 ): ¯ ¯ ¯ P K ¡ 1 k =0 j ^ ® i ( ! k ) j 2 ¡ P K ¡ 1 k =0 j ^ ® i ¡ 1 ( ! k ) j 2 ¯ ¯ ¯ P K ¡ 1 k =0 j ^ ® i ¡ 1 ( ! k ) j 2 · ²; (3.51) where w e ev aluate ^ ® ( ! ) on a K -p oin t DFT grid: ! k = 2 ¼ k =K for k = 0 ; :::; K ¡ 1. 3.6 MAPES Compared with GAPES As explained ab o v e, MAPES is deriv ed from a surrogate ML form ulation of the APES algorithm; on the other hand, GAPES is deriv ed from a least squares (LS) form ulation of APES [69]. In the complete-data case, these t w o approac hes are equiv alen t in the sense that from either of them w e can deriv e the same fulldata APES sp ectral estimator. So at ¯rst, it migh t lo ok coun terin tuitiv e that these t w o algorithms (MAPES and GAPES) will p erform di®eren tly for the missing-data problem (see the n umerical results in Section 5.7). Here w e will giv e a brief discussion ab out this issue. The di®erence b et w een MAPES and GAPES concerns the w a y they estimate ¹ when some data samples are missing. (See App endix 3.9 for a detailed deriv ation

PAGE 43

31 of GAPES.) Although MAPES-EM estimates eac h missing sample separately for eac h frequency ! k (and for eac h data snapshot ¹ y l in MAPES-EM2) while GAPES estimates eac h missing sample b y considering all K frequencies together, the real di®erence b et w een them concerns the di®eren t criteria used in (3.68) and (3.3) for the estimation of ¹ : GAPES estimates the missing sample ¹ based on a least squares (LS) ¯tting of the ¯ltered data, h H ( ! k ) ¹ y l . On the other hand, MAPES estimates the missing sub-sample ¹ directly from f ¹ y l g based on an \ML" ¯tting criterion. Due to the fact that the LS form ulation of APES fo cuses on the output of the ¯lter h ( ! k ) (whic h is supp osed to suppress an y other frequency comp onen ts except ! k ), the GAPES algorithm is sensitiv e to the errors in h ( ! k ) when it tries to estimate the missing data. This is the reason wh y GAPES p erforms w ell in the gapp ed-data case, since there a go o d estimate of h ( ! k ) can b e calculated during the initialization step. When the missing samples o ccur in an arbitrary pattern, ho w ev er, the p erformance of GAPES degrades. Y et the MAPES-EM do es not su®er from suc h a degradation. 3.7 Numerical Examples W e presen t detailed results of a few n umerical examples to demonstrate the p erformance of the MAPES-EM algorithms for missing-data sp ectral estimation. W e compare MAPES-EM with WFFT and GAPES. A T a ylor windo w with order 5 and sidelob e lev el -35 dB is used for WFFT. W e c ho ose K = 32 N to ha v e a ¯ne grid of discrete frequencies. W e calculate the corresp onding WFFT sp ectrum via zero-padded FFT. The so-obtained WFFT sp ectrum is used as the initial sp ectral estimate for the MAPES-EM and GAPES algorithms. The initial estimate of Q ( ! ) for MAPES-EM has b een discussed b efore, and the initial estimate of h ( ! ) for GAPES is calculated from (2.12), where the missing samples are set to zero. W e stop the MAPES-EM and the GAPES algorithms using the same stopping criterion in (3.51) with ² b eing selected as 10 ¡ 3 and 10 ¡ 2 , resp ectiv ely . The reason w e c ho ose a larger ² for GAPES is that it con v erges relativ ely slo wly for the general missing-data problem and its

PAGE 44

32 sp ectral estimate w ould not impro v e m uc h if w e used an ² < 10 ¡ 2 . All the adaptiv e ¯ltering algorithms considered (i.e. APES, GAPES, and MAPES-EM) use a ¯lter length of M = N = 2 for ac hieving high resolution. The true sp ectrum of the sim ulated signal is sho wn in Figure 3.1(a), where w e ha v e four sp ectral lines lo cated at f 1 = 0 : 05 Hz, f 2 = 0 : 065 Hz, f 3 = 0 : 26 Hz, and f 4 = 0 : 28 Hz with complex amplitudes ® 1 = ® 2 = ® 3 = 1 and ® 4 = 0 : 5. Besides these sp ectral lines, Figure 3.1(a) also sho ws a con tin uous sp ectral comp onen t cen tered at 0.18 Hz with a width b = 0 : 015 Hz and a constan t mo dulus of 0.25. The data sequence has N = 128 samples among whic h 51 (40%) samples are missing; the lo cations of the missing samples are c hosen arbitrarily . The data is corrupted b y a zero-mean circularly symmetric complex white Gaussian noise with v ariance ¾ 2 n = 0 : 01. In Figure 3.1(b), the APES algorithm is applied to the complete data and the resulting sp ectrum is sho wn. The APES sp ectrum will b e used later as a reference for comparison purp oses. The WFFT sp ectrum for the incomplete data is sho wn in Figure 3.1(c), where the artifacts due to the missing data are readily observ ed. As exp ected, the WFFT sp ectrum has p o or resolution and high sidelob es and it underestimates the true sp ectrum. Note that the WFFT sp ectrum will b e used as the initial estimate for the GAPES and MAPES algorithms. Figure 3.1(d) sho ws the GAPES sp ectrum. GAPES also underestimates the sin usoidal comp onen ts and giv es some artifacts. Apparen tly , due to the p o or initial estimate of h ( ! k ) for the incomplete data, GAPES con v erges to one of the lo cal minima of the cost function in (3.68). Figures. 3.1(e) and 3.1(f ) sho w the MAPES-EM1 and MAPES-EM2 sp ectral estimates. Both MAPES algorithms p erform quite w ell and their sp ectral estimates are similar to the high-resolution APES sp ectrum in Figure 3.1(b). The MAPES-EM1 and MAPES-EM2 sp ectral estimates at di®eren t iterations are plotted in Figures. 3.2(a) and 3.2(b), resp ectiv ely . Both algorithms con v erge

PAGE 45

33 quic kly with MAPES-EM1 con v erging after 10 iterations while MAPES-EM2 after only 6. The data restoration p erformance of MAPES-EM is sho wn in Figure 3.3. The missing samples are estimated using the a v eraging approac h w e in tro duced previously . Figures. 3.3(a) and 3.3(b) displa y the real and imaginary parts of the in terp olated data, resp ectiv ely , obtained via MAPES-EM1. Figures. 3.3(c) and 3.3(d) sho w the corresp onding results for MAPES-EM2. The lo cations of the missing samples are also indicated in Figure 3.3. The missing samples estimated via the MAPES-EM algorithms are quite accurate. More detailed results for MAPES-EM2 are sho wn in Fig 3.4. (Those for MAPES-EM1 are similar.) F or a clear visualization, only the estimates of the ¯rst three missing samples are sho wn in Figure 3.4. The real and imaginary parts of the estimated samples as a function of frequency are plotted in Figures. 3.4(a) and 3.4(b), resp ectiv ely . All estimates are close to the corresp onding true v alues whic h are also indicated in Figure 3.4. It is in teresting to note that larger v ariations o ccur at frequencies where strong signal comp onen ts are presen t. The results displa y ed so far w ere for one randomly pic k ed realization of the data. Using 100 Mon te Carlo sim ulations (v arying the realizations of the noise, the initial phases of the di®eren t sp ectral comp onen ts, and the missing-data patterns), w e obtain the ro ot mean-squared errors (RMSE's) of the magnitude and phase estimates of the four sp ectral lines at their true frequency lo cations. These RMSE's for WFFT, GAPES, and MAPES-EM are listed in T ables 3.1 and 3.2. Based on this limited set of Mon te-Carlo sim ulations, w e can see that the t w o MAPES-EM algorithms p erform similarly , and that they are m uc h more accurate than WFFT and GAPES. A similar b eha vior has b een observ ed in sev eral other n umerical exp erimen ts. Next, w e increase the n um b er of missing samples to 77 (60% of the original data). The results of WFFT, WFFT GAPES, MAPES-EM1, and MAPES-EM2 are

PAGE 46

34 T able 3.1: RMSE's of the magnitude estimates obtained via the WFFT, GAPES, and MAPES-EM sp ectral estimators. WFFT GAPES MAPES-EM1 MAPES-EM2 Signal 1 0.420 0.175 0.010 0.011 Signal 2 0.417 0.205 0.010 0.010 Signal 3 0.419 0.169 0.010 0.009 Signal 4 0.205 0.164 0.009 0.011 T able 3.2: RMSE's of the phase (radian) estimates obtained via the WFFT, GAPES, and MAPES-EM sp ectral estimators. WFFT GAPES MAPES-EM1 MAPES-EM2 Signal 1 0.077 0.042 0.021 0.021 Signal 2 0.059 0.038 0.024 0.025 Signal 3 0.099 0.037 0.012 0.013 Signal 4 0.133 0.029 0.022 0.022 sho wn in Figures. 3.5(a), 3.5(b), 3.5(c), and 3.5(d), resp ectiv ely . The signal amplitudes in the WFFT sp ectrum are lo w presumably due to the small (40%) amoun t of a v ailable data samples. The artifacts are so high that w e can hardly iden tify the signal comp onen ts. GAPES also p erforms p o orly (see Figure 3.5(b)), as exp ected. The MAPES-EM algorithms pro vide excellen t sp ectral estimates with relativ ely minor artifacts. In our next exp erimen t, w e k eep the 40% data missing rate but increase the noise v ariance to ¾ 2 n = 0 : 1 (10 dB higher than in the previous exp erimen ts). The corresp onding mo duli of the sp ectral estimates of complete-data WFFT, APES, missing-data WFFT, GAPES, MAPES-EM1, and MAPES-EM2 are plotted in Figures. 3.6(a), 3.6(b), 3.6(c), 3.6(d), 3.6(e), and 3.6(f ), resp ectiv ely . Again, the p erformance of the MAPES-EM algorithms is excellen t. In our last exp erimen t, w e plot the RMSE's of the MAPES-EM1 estimates as a function of the missing sample rate in Figure 3.7. F or clear illustration, only the

PAGE 47

35 RMSE's of the estimates of ¯rst sp ectral line lo cated at f 1 = 0 : 05 Hz are plotted. Eac h result is based on 100 Mon te Carlo sim ulations (v arying the realizations of the noise, the initial phases of the di®eren t sp ectral comp onen ts, and the missing-data patterns). Here the signal-to-noise ratio (SNR) is de¯ned as SNR 1 = 10 log 10 j ® 1 j 2 ¾ 2 n (dB) ; (3.52) with ® 1 b eing the complex amplitude of the ¯rst sin usoid. Clearly , for eac h ¯xed SNR, the RMSE's increase as the n um b er of missing samples increases. Also as exp ected, the RMSE's decrease when the SNR increases. Similar results can b e observ ed via MAPES-EM2. 3.8 Conclusions W e ha v e considered nonparametric complex sp ectral estimation for a data sequence with missing samples o ccurring in arbitrary patterns. Tw o MAPES-EM algorithms, namely MAPES-EM1 and MAPES-EM2, ha v e b een deriv ed b y form ulating an ML based problem whic h is solv ed iterativ ely via t w o EM algorithms. F urthermore, an a v eraging approac h has also b een suggested for missing-sample estimation via MAPES-EM. W e ha v e compared MAPES-EM with the previously dev elop ed GAPES algorithm and demonstrated the adv an tages of the prop osed algorithms. Numerical sim ulations ha v e b een pro vided to demonstrate the excellen t p erformance of the t w o MAPES-EM algorithms for missing-data sp ectral estimation and missing data restoration. The t w o MAPES-EM algorithms ha v e p erformed quite similarly in sim ulations, but MAPES-EM2 is computationally more app ealing than MAPESEM1. 3.9 App endix: GAPES 3.9.1 In tro duction One sp ecial case of the missing-data problem is called gapp ed-data, where the measuremen ts during certain p erio ds are not v alid due to man y reasons suc h

PAGE 48

36 0 0.1 0.2 0.3 0.4 0.5 0 0.2 0.4 0.6 0.8 1 1.2 Frequency (Hz)Modulus of Complex Amplitude 0 0.1 0.2 0.3 0.4 0.5 0 0.2 0.4 0.6 0.8 1 1.2 Frequency (Hz)Modulus of Complex Amplitude (a) (b) 0 0.1 0.2 0.3 0.4 0.5 0 0.2 0.4 0.6 0.8 1 1.2 Frequency (Hz)Modulus of Complex Amplitude 0 0.1 0.2 0.3 0.4 0.5 0 0.2 0.4 0.6 0.8 1 1.2 Frequency (Hz)Modulus of Complex Amplitude (c) (d) 0 0.1 0.2 0.3 0.4 0.5 0 0.2 0.4 0.6 0.8 1 1.2 Frequency (Hz)Modulus of Complex Amplitude 0 0.1 0.2 0.3 0.4 0.5 0 0.2 0.4 0.6 0.8 1 1.2 Frequency (Hz)Modulus of Complex Amplitude (e) (f ) Figure 3.1: Mo dulus of the missing-data sp ectral estimates ( N = 128, ¾ 2 n = 0 : 01, 51 (40%) missing samples): (a) T rue sp ectrum, (b) Complete-data APES, (c) WFFT, (d) GAPES with M = 64 and ² = 10 ¡ 2 , (e) MAPES-EM1 with M = 64 and ² = 10 ¡ 3 , and (f ) MAPES-EM2 with M = 64 and ² = 10 ¡ 3 .

PAGE 49

37 0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 Modulus i=0 0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 Modulus i=1 0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 Modulus i=2 0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 Modulus i=3 0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 Modulus i=4 0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 Modulus i=5 Frequency (Hz) 0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 Modulus i=0 0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 Modulus i=1 0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 Modulus i=2 0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 Modulus i=3 0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 Modulus i=4 0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 Modulus i=5 Frequency (Hz) (a) (b) Figure 3.2: Mo dulus of the missing-data sp ectral estimates obtained via the MAPESEM algorithms at di®eren t iterations ( N = 128, ¾ 2 n = 0 : 01, 51 (40%) missing samples): (a) MAPES-EM1, and (b) MAPES-EM2.

PAGE 50

38 20 40 60 80 100 120 -4 -2 0 2 4 nReal Part of Signal True Data Interpolated Data Missing Data Locations (a) 20 40 60 80 100 120 -4 -2 0 2 4 nImaginary Part of Signal True Data Interpolated Data Missing Data Locations (b) 20 40 60 80 100 120 -4 -2 0 2 4 nReal Part of Signal True Data Interpolated Data Missing Data Locations (c) 20 40 60 80 100 120 -4 -2 0 2 4 nImaginary Part of Signal True Data Interpolated Data Missing Data Locations (d) Figure 3.3: In terp olation of the missing samples ( N = 128, ¾ 2 n = 0 : 01, 51 (40%) missing samples): (a) Real part of the data in terp olated via MAPES-EM1, (b) Imaginary part of the data in terp olated via MAPES-EM1, (c) Real part of the data in terp olated via MAPES-EM2, and (d) Imaginary part of the data in terp olated via MAPES-EM2.

PAGE 51

39 0 0.1 0.2 0.3 0.4 0.5 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 Frequency (Hz)Real Part of Missing Sample Estimates 1st Missing Sample 2nd Missing Sample 3rd Missing Sample True Value (1st) True Value (2nd) True Value (3rd) (a) 0 0.1 0.2 0.3 0.4 0.5 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 Frequency (Hz)Imaginary Part of Missing Sample Estimates 1st Missing Sample 2nd Missing Sample 3rd Missing Sample True Value (1st) True Value (2nd) True Value (3rd) (b) Figure 3.4: Estimation of the ¯rst three missing samples vs. frequency ( N = 128, ¾ 2 n = 0 : 01, 51 (40%) samples are missing, v ertical dotted lines indicate the true frequency lo cations of the sp ectral comp onen ts with the closely spaced lines indicating the con tin uous sp ectral comp onen t): (a) Real part and (b) Imaginary part of the missing samples estimated via MAPES-EM2.

PAGE 52

40 0 0.1 0.2 0.3 0.4 0.5 0 0.2 0.4 0.6 0.8 1 1.2 Frequency (Hz)Modulus of Complex Amplitude 0 0.1 0.2 0.3 0.4 0.5 0 0.2 0.4 0.6 0.8 1 1.2 Frequency (Hz)Modulus of Complex Amplitude (a) (b) 0 0.1 0.2 0.3 0.4 0.5 0 0.2 0.4 0.6 0.8 1 1.2 Frequency (Hz)Modulus of Complex Amplitude 0 0.1 0.2 0.3 0.4 0.5 0 0.2 0.4 0.6 0.8 1 1.2 Frequency (Hz)Modulus of Complex Amplitude (c) (d) Figure 3.5: Mo dulus of the missing-data sp ectral estimates ( N = 128, ¾ 2 n = 0 : 01, 77 (60%) missing samples) obtained via: (a) WFFT, (b) GAPES with M = 64 and ² = 10 ¡ 2 , (c) MAPES-EM1 with M = 64 and ² = 10 ¡ 3 , and (d) MAPES-EM2 with M = 64 and ² = 10 ¡ 3 .

PAGE 53

41 0 0.1 0.2 0.3 0.4 0.5 0 0.2 0.4 0.6 0.8 1 1.2 Frequency (Hz)Modulus of Complex Amplitude 0 0.1 0.2 0.3 0.4 0.5 0 0.2 0.4 0.6 0.8 1 1.2 Frequency (Hz)Modulus of Complex Amplitude (a) (b) 0 0.1 0.2 0.3 0.4 0.5 0 0.2 0.4 0.6 0.8 1 1.2 Frequency (Hz)Modulus of Complex Amplitude 0 0.1 0.2 0.3 0.4 0.5 0 0.2 0.4 0.6 0.8 1 1.2 Frequency (Hz)Modulus of Complex Amplitude (c) (d) 0 0.1 0.2 0.3 0.4 0.5 0 0.2 0.4 0.6 0.8 1 1.2 Frequency (Hz)Modulus of Complex Amplitude 0 0.1 0.2 0.3 0.4 0.5 0 0.2 0.4 0.6 0.8 1 1.2 Frequency (Hz)Modulus of Complex Amplitude (e) (f ) Figure 3.6: Mo dulus of the missing-data sp ectral estimates ( N = 128, ¾ 2 n = 0 : 1, 51 (40%) missing samples) obtained via: (a) Complete-data WFFT, (b) Complete-data APES, (c) WFFT, (d) GAPES with M = 64 and ² = 10 ¡ 2 , (e) MAPES-EM1 with M = 64 and ² = 10 ¡ 3 , and (f ) MAPES-EM2 with M = 64 and ² = 10 ¡ 3 .

PAGE 54

42 0.1 0.2 0.3 0.4 0.5 0.6 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 missing sample rateRMSE SNR=0 SNR=5 SNR=10 SNR=15 SNR=20 (a) 0.1 0.2 0.3 0.4 0.5 0.6 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 missing sample rateRMSE SNR=0 SNR=5 SNR=10 SNR=15 SNR=20 (b) Figure 3.7: RMSE's of the estimates of the ¯rst sp ectral line lo cated at f 1 = 0 : 05 Hz obtained via MAPES-EM1: (a) amplitude, and (b) phase (radian).

PAGE 55

43 as in terference or jamming. The di®erence b et w een the gapp ed-data problem and the general missing-data problem, where the missing samples can o ccur at arbitrary places among the complete data set, is that for the gapp ed-data case, there exists group(s) of a v ailable data samples where within eac h group, there are no missing samples. Suc h scenarios exist in astronomical or radar applications where large segmen ts of data are a v ailable in spite of the fact that the data b et w een these segmen ts are missing. F or example, in radar signal pro cessing, the problem of com bining several sets of measuremen ts made at di®eren t azim uth angle lo cations can b e p osed as a problem of sp ectral estimation from gapp ed data [30 , 29 ]. Similar problems arise in data fusion via ultra wideband coheren t pro cessing [12 ]. In astronom y , data are often a v ailable as groups of samples with rather long in terv als during whic h no measuremen ts can b e tak en [44 , 62 , 63, 53 , 59, 51 ]. The gapp ed-data APES (GAPES) considers using the APES ¯lter (as in troduced in Chapter 2) for the sp ectral estimation of gapp ed-data. Sp eci¯cally , the GAPES algorithm consists of t w o steps: 1) estimating the adaptiv e ¯lter and the corresp onding sp ectrum via APES, and 2) ¯lling in the gaps via a least squares (LS) ¯t. 3.9.2 GAPES Assume that some segmen ts of the 1-D data sequence f y n g N ¡ 1 n =0 are una v ailable. Let y 4 = [ y 1 y 2 ¢ ¢ ¢ y N ¡ 1 ] T 4 = [ y T 1 y T 2 ¢ ¢ ¢ y T P ] T (3.53) b e the complete data v ector, where y 1 ; :::; y P are sub v ectors of y , whose lengths are N 1 ; :::; N P , resp ectiv ely , with N 1 + N 2 + ¢ ¢ ¢ + N P = N . A gapp ed-data v ector ° is

PAGE 56

44 formed b y assuming y p , for p = 1 ; 3 ; :::; P ( P is alw a ys an o dd n um b er), are a v ailable: ° 4 = [ y T 1 y T 3 ¢ ¢ ¢ y T P ] T : (3.54) Similarly , ¹ 4 = [ y T 2 y T 4 ¢ ¢ ¢ y T P ¡ 1 ] T (3.55) denotes all the missing samples. Then ° and ¹ ha v e dimensions g £ 1 and ( N ¡ g ) £ 1, resp ectiv ely , where g = N 1 + N 3 + ¢ ¢ ¢ + N P is the total n um b er of a v ailable samples. 3.9.2.1 Initial Estimates via APES W e obtain the initial APES estimates of h ( ! ) and ® ( ! ) from the a v ailable data ° as follo ws. Cho ose an initial ¯lter length M 0 suc h that an initial full-rank co v ariance matrix ^ R can b e built with the ¯lter length M 0 using the a v ailable data segmen ts only . This indicates X p 2f 1 ; 3 ;:::;P g max(0 ; N p ¡ M 0 + 1) > M 0 : (3.56) Let L p = N p ¡ M 0 + 1 and let J b e the subset of f 1 ; 3 ; :::; P g for whic h L p > 0. Then the ¯lter-bank h ( ! ) is calculated from (2.11) and (2.12) b y using the follo wing rede¯nitions: ^ R = 1 P p 2J L p X p 2J N 1 + ¢¢¢ + N p ¡ 1 + L p ¡ 1 X l = N 1 + ¢¢¢ + N p ¡ 1 ¹ y l ¹ y H l ; (3.57) ¹ g ( ! ) = 1 P p 2J L p X p 2J N 1 + ¢¢¢ + N p ¡ 1 + L p ¡ 1 X l = N 1 + ¢¢¢ + N p ¡ 1 ¹ y l e ¡ j ! l : (3.58) Note that the data snapshots used ab o v e ha v e a size of M 0 £ 1 whose elemen ts are from ° only , and hence they do not con tain an y missing samples. Corresp ondingly , the ^ R and ¹ g ( ! ) estimated ab o v e ha v e sizes of M 0 £ M 0 and M 0 £ 1, resp ectiv ely . Next, the ¯lter-bank h ( ! ) is applied to the a v ailable data ° and the LS estimate of ® ( ! ) from the ¯lter output is calculated b y using (2.16) where ¹ g ( ! ) is

PAGE 57

45 replaced b y (3.58). Note that in the ab o v e ¯ltering pro cess, only a v ailable samples are passed through the ¯lter. The initial LS estimate of ® ( ! ) is based on these so-obtained ¯lter outputs only . 3.9.2.2 Data In terp olation No w w e consider the estimation of ¹ based on the initial sp ectral estimates ^ ® ( ! ) and ^ h ( ! ) obtained as outlined ab o v e. Under the assumption that the missing data ha v e the same sp ectral con ten t as the a v ailable data, w e can determine ^ ° under the condition that the output of the ¯lter ^ h ( ! ) fed with the complete data sequence made from ° and ^ ¹ is as close as p ossible (in the LS sense) to ^ ® ( ! ) e ¡ j ! l for l = 0 ; :::; L ¡ 1. Since usually w e ev aluate ^ ® ( ! ) on a K -p oin t DFT grid: ! k = 2 ¼ k =K for k = 0 ; :::; K ¡ 1 (usually w e ha v e K > N ), w e obtain ¹ as the solution to the follo wing LS problem: min ¹ K ¡ 1 X k =0 L ¡ 1 X l =0 ¯ ¯ ¯ ^ h H ( ! k ) ¹ y l ¡ ^ ® ( ! k ) e j ! k l ¯ ¯ ¯ 2 (3.59) Note that b y estimating ¹ this w a y , w e remain in the LS ¯tting framew ork of APES. The quadratic minimization problem (3.59) can b e readily solv ed. Let H ( ! k ) = 2 6 6 6 4 ^ h ¤ 0 ¢ ¢ ¢ ^ h ¤ M 0 ¡ 1 0 0 0 0 ^ h ¤ 0 ¢ ¢ ¢ ^ h ¤ M 0 ¡ 1 0 0 . . . . . . . . . 0 0 0 ^ h ¤ 0 ¢ ¢ ¢ ^ h ¤ M 0 ¡ 1 3 7 7 7 5 = 2 6 6 6 4 ^ h H ( ! k ) ^ h H ( ! k ) . . . ^ h H ( ! k ) 3 7 7 7 5 2 C L £ N (3.60) and ´ ( ! k ) = ^ ® ( ! k ) 2 6 6 6 4 1 e j ! k . . . e j ! k ( L ¡ 1) 3 7 7 7 5 2 C L £ 1 : (3.61) Using this notation w e can write the ob jectiv e function in (3.59) as K ¡ 1 X k =0 ° ° ° ° ° ° ° H ( ! k ) 2 6 4 y 0 . . . y N ¡ 1 3 7 5 ¡ ´ ( ! k ) ° ° ° ° ° ° ° 2 (3.62)

PAGE 58

46 De¯ne the L £ g and L £ ( N ¡ g ) matrices A ( ! k ) and B ( ! k ) from H ( ! k ) via the follo wing equalit y: H ( ! k ) 2 6 4 y 0 . . . y N ¡ 1 3 7 5 = A ( ! k ) ° + B ( ! k ) ¹ : (3.63) Also, let d ( ! k ) = ´ ( ! k ) ¡ A ( ! k ) ° : (3.64) With this notation the ob jectiv e function (3.62) b ecomes K ¡ 1 X k =0 k B ( ! k ) ¹ ¡ d ( ! k ) k 2 ; (3.65) whose minimizer with resp ect to ¹ is readily found to b e ^ ¹ = Ã K ¡ 1 X k =0 B H ( ! k ) B ( ! k ) ! ¡ 1 Ã K ¡ 1 X k =0 B H ( ! k ) d ( ! k ) ! : (3.66) 3.9.2.3 Summary of GAPES Once an estimate ^ ¹ has b ecome a v ailable, the next logical step should consist of re-estimating the sp ectrum and the ¯lter-bank, b y applying APES to the data sequence made from ° and ^ ¹ . According to the discussion around (2.4), this en tails the minimization with resp ect to h ( ! k ) and ® ( ! k ) of the function K ¡ 1 X k =0 L ¡ 1 X l =0 ¯ ¯ h H ( ! k ) ^ ¹ y l ¡ ® ( ! k ) e j ! k l ¯ ¯ 2 (3.67) sub ject to h H ( ! k ) a ( ! k ) = 1, where ^ ¹ y l is made from ° and ^ ¹ . Eviden tly , the minimization of (3.67) with resp ect to f h ( ! k ) ; ® ( ! k ) g K ¡ 1 k =0 can b e decoupled in to K minimization problems of the form of (2.4), y et w e prefer to write the criterion function as in (3.67) to mak e the connection with (3.59). In e®ect, comparing (3.59) and (3.67) clearly sho ws that the alternating estimation of f ® ( ! k ) ; h ( ! k ) g and ¹ outlined ab o v e can b e recognized as a cyclic optimization approac h for solving the follo wing minimization problem: min ¹ ; f ® ( ! k ) ; h ( ! k ) g K ¡ 1 X k =0 L ¡ 1 X l =0 ¯ ¯ h H ( ! k ) ¹ y l ¡ ® ( ! k ) e j ! k l ¯ ¯ 2 sub ject to h H ( ! k ) a ( ! k ) = 1 ; (3.68)

PAGE 59

47 A step-b y-step summary of GAPES follo ws. Step 0: Obtain an initial estimate of f ® ( ! k ) ; h ( ! k ) g . Step 1: Use the most recen t estimate of f ® ( ! k ) ; h ( ! k ) g in (3.68) to estimate ¹ b y minimizing the so-obtained cost function, whose solution is giv en b y (3.66). Step 2: Use the latest estimate of ¹ to ¯ll in the missing data samples, and estimate f ® ( ! k ) ; h ( ! k ) g K ¡ 1 k =0 b y minimizing the cost function in (3.68) based on the in terp olated data. (This step is equiv alen t to applying APES to the complete data.) Step 3: Rep eat Steps 1 to 2 un til practical con v ergence. The practical con v ergence can b e decided when the relativ e c hange of the cost function in (3.68) corresp onding to the curren t and previous estimates is smaller than a preassigned threshold (e.g., ² = 10 ¡ 3 ). After con v ergence, w e ha v e a ¯nal sp ectral estimate f ^ ® ( ! k ) g K ¡ 1 k =0 . If desired, w e can use the ¯nal in terp olated data sequence to compute the APES sp ectrum on a grid ev en ¯ner than the one used in the aforemen tioned minimization pro cedure. Note that usually the selected initial ¯lter length satis¯es M 0 < M due to the missing data samples, so there are man y practical c hoices to increase the ¯lter length after initialization, whic h include, for example, increasing the ¯lter length after eac h iteration un til it reac hes M . F or simplicit y , w e c ho ose to use ¯lter length M righ t after the initialization step.

PAGE 60

CHAPTER 4 TW O-DIMENSIONAL APES F OR COMPLETE D A T A SPECTRAL ESTIMA TION 4.1 In tro duction In the previous c hapters, 1-D APES and MAPES algorithms ha v e b een in troduced for complete-data and missing-data sp ectral estimation. In this c hapter, w e pro vide the 2-D extension of the APES algorithm. The deriv ations of the 2-D APES sp ectral estimator based on the adaptiv e ¯ltering and the maxim um lik eliho o d ¯tting approac hes are presen ted in Sections 4.2 and 4.3, resp ectiv ely . 4.2 2-D APES Filter Consider the problem of estimating the amplitude sp ectrum of a complexv alued uniformly sampled 2-D discrete-time signal f y n 1 ;n 2 g N 1 ¡ 1 ;N 2 ¡ 1 n 1 =0 ;n 2 =0 , where the data matrix has dimension N 1 £ N 2 . F or a 2-D frequency ( ! 1 ; ! 2 ) of in terest, the signal y n 1 ;n 2 is describ ed as y n 1 ;n 2 = ® ( ! 1 ; ! 2 ) e j ( ! 1 n 1 + ! 2 n 2 ) + e n 1 ;n 2 ( ! 1 ; ! 2 ) ; (4.1) n 1 = 0 ; :::; N 1 ¡ 1 ; n 2 = 0 ; :::; N 2 ¡ 1 ; ! 1 ; ! 2 2 [0 ; 2 ¼ ) ; where ® ( ! 1 ; ! 2 ) denotes the complex amplitude of the 2-D sin usoidal comp onen t at frequency ( ! 1 ; ! 2 ), and e n 1 ;n 2 ( ! 1 ; ! 2 ) denotes the residual matrix (assumed zeromean) whic h includes the unmo deled noise and in terference from frequencies other than ( ! 1 ; ! 2 ). The 2-D APES algorithm deriv ed b elo w estimates ® ( ! 1 ; ! 2 ) from f y n 1 ;n 2 g for an y giv en frequency pair ( ! 1 ; ! 2 ). 48

PAGE 61

49 Let Y b e an N 1 £ N 2 data matrix Y 4 = 2 6 6 6 4 y 0 ; 0 y 0 ; 1 ¢ ¢ ¢ y 0 ;N 2 ¡ 1 y 1 ; 0 y 1 ; 1 ¢ ¢ ¢ y 1 ;N 2 ¡ 1 . . . . . . . . . . . . y N 1 ¡ 1 ; 0 y N 1 ¡ 1 ; 1 ¢ ¢ ¢ y N 1 ¡ 1 ;N 2 ¡ 1 3 7 7 7 5 ; (4.2) and let H ( ! 1 ; ! 2 ) b e an M 1 £ M 2 matrix that con tains the co e±cien ts of a 2-D FIR ¯lter H ( ! 1 ; ! 2 ) 4 = 2 6 6 6 4 h 0 ; 0 ( ! 1 ; ! 2 ) h 0 ; 1 ( ! 1 ; ! 2 ) ¢ ¢ ¢ h 0 ;M 2 ¡ 1 ( ! 1 ; ! 2 ) h 1 ; 0 ( ! 1 ; ! 2 ) h 1 ; 1 ( ! 1 ; ! 2 ) ¢ ¢ ¢ h 1 ;M 2 ¡ 1 ( ! 1 ; ! 2 ) . . . . . . . . . . . . h M 1 ¡ 1 ; 0 ( ! 1 ; ! 2 ) h M 1 ¡ 1 ; 1 ( ! 1 ; ! 2 ) ¢ ¢ ¢ h M 1 ¡ 1 ;M 2 ¡ 1 ( ! 1 ; ! 2 ) 3 7 7 7 5 : (4.3) Let L 1 4 = N 1 ¡ M 1 + 1 and L 2 4 = N 2 ¡ M 2 + 1. Then w e denote b y ¹ X = H ( ! 1 ; ! 2 ) ? Y (4.4) the follo wing L 1 £ L 2 output data matrix obtained b y ¯ltering Y through the ¯lter determined b y H ( ! 1 ; ! 2 ) ¹ x l 1 ;l 2 = M 1 ¡ 1 X m 1 =0 M 2 ¡ 1 X m 2 =0 h ¤ m 1 ;m 2 ( ! 1 ; ! 2 ) y l 1 + m 1 ;l 2 + m 2 = v ec H ( H ( ! 1 ; ! 2 )) ¹ y l 1 ;l 2 ; (4.5) where v ec( ¢ ) denotes the op eration of stac king the columns of a matrix on to eac h other. In (4.5), ¹ y l 1 ;l 2 is de¯ned b y ¹ y l 1 ;l 2 4 = v ec( ¹ Y l 1 ;l 2 ) 4 = v ec 0 B B B @ 2 6 6 6 4 y l 1 ;l 2 y l 1 ;l 2 +1 ¢ ¢ ¢ y l 1 ;l 2 + M 2 ¡ 1 y l 1 +1 ;l 2 y l 1 +1 ;l 2 +1 ¢ ¢ ¢ y l 1 +1 ;l 2 + M 2 ¡ 1 . . . . . . . . . . . . y l 1 + M 1 ¡ 1 ;l 2 y l 1 + M 1 ¡ 1 ;l 2 +1 ¢ ¢ ¢ y l 1 + M 1 ¡ 1 ;l 2 + M 2 ¡ 1 3 7 7 7 5 1 C C C A : (4.6) The APES sp ectrum estimate ^ ® ( ! 1 ; ! 2 ) and the ¯lter co e±cien t matrix ^ H ( ! 1 ; ! 2 ) are the minimizers of the follo wing LS criterion: min ® ( ! 1 ;! 2 ) ; H ( ! 1 ;! 2 ) L 1 ¡ 1 X l 1 =0 L 2 ¡ 1 X l 2 =0 ¯ ¯ ¹ x l 1 ;l 2 ¡ ® ( ! 1 ; ! 2 ) e j ( ! 1 l 1 + ! 2 l 2 ) ¯ ¯ 2 sub ject to v ec H ( H ( ! 1 ; ! 2 )) a M 1 ;M 2 ( ! 1 ; ! 2 ) = 1 : (4.7)

PAGE 62

50 Here a M 1 ;M 2 ( ! 1 ; ! 2 ) is an M 1 M 2 £ 1 v ector giv en b y a M 1 ;M 2 ( ! 1 ; ! 2 ) 4 = a M 2 ( ! 2 ) ­ a M 1 ( ! 1 ) ; (4.8) where ­ denotes the Kronec k er matrix pro duct, and a M k ( ! k ) 4 = [1 e j ! k ¢ ¢ ¢ e j ( M k ¡ 1) ! k ] T ; k = 1 ; 2 : (4.9) Substituting ¹ X in to (4.7), w e ha v e the follo wing design ob jectiv e for 2-D APES: min ® ( ! 1 ;! 2 ) ; H ( ! 1 ;! 2 ) k v ec ( H ( ! 1 ; ! 2 ) ? Y ) ¡ ® ( ! 1 ; ! 2 ) a L 1 ;L 2 ( ! 1 ; ! 2 ) k 2 sub ject to v ec H ( H ( ! 1 ; ! 2 )) a M 1 ;M 2 ( ! 1 ; ! 2 ) = 1 ; (4.10) where a L 1 ;L 2 ( ! 1 ; ! 2 ) is de¯ned similar to that of a M 1 ;M 2 ( ! 1 ; ! 2 ). The solution to (4.10) can b e readily deriv ed. De¯ne ^ R = 1 L 1 L 2 L 1 ¡ 1 X l 1 =0 L 2 ¡ 1 X l 2 =0 ¹ y l 1 ;l 2 ¹ y H l 1 ;l 2 (4.11) and let ¹ g ( ! 1 ; ! 2 ) denote the normalized 2-D F ourier transform of ¹ y l 1 ;l 2 : ¹ g ( ! 1 ; ! 2 ) = 1 L 1 L 2 L 1 ¡ 1 X l 1 =0 L 2 ¡ 1 X l 2 =0 ¹ y l 1 ;l 2 e ¡ j ( ! 1 l 1 + ! 2 l 2 ) : (4.12) The ¯lter H ( ! 1 ; ! 2 ) that minimizes (4.10) is giv en b y v ec ( ^ H ( ! 1 ; ! 2 )) = ^ S ¡ 1 ( ! 1 ; ! 2 ) a M 1 ;M 2 ( ! 1 ; ! 2 ) a H M 1 ;M 2 ( ! 1 ; ! 2 ) ^ S ¡ 1 ( ! 1 ; ! 2 ) a M 1 ;M 2 ( ! 1 ; ! 2 ) (4.13) and the APES sp ectrum is giv en b y ^ ® ( ! 1 ; ! 2 ) = a H M 1 ;M 2 ( ! 1 ; ! 2 ) ^ S ¡ 1 ( ! 1 ; ! 2 ) ¹ g ( ! 1 ; ! 2 ) a H M 1 ;M 2 ( ! 1 ; ! 2 ) ^ S ¡ 1 ( ! 1 ; ! 2 ) a M 1 ;M 2 ( ! 1 ; ! 2 ) ; (4.14) where ^ S ( ! 1 ; ! 2 ) 4 = ^ R ¡ ¹ g ( ! 1 ; ! 2 ) ¹ g H ( ! 1 ; ! 2 ) : (4.15)

PAGE 63

51 4.3 2-D ML-Based APES W e pro vide herein the 2-D extension of APES, devised via a \ML ¯tting" based approac h, for complete-data sp ectral estimation. Consider the 2-D problem in tro duced in Section 4.2. P artition the N 1 £ N 2 data matrix Y in to L 1 L 2 o v erlapping submatrices of size M 1 £ M 2 : ¹ Y l 1 ;l 2 = 2 6 6 6 4 y l 1 ;l 2 y l 1 ;l 2 +1 ¢ ¢ ¢ y l 1 ;l 2 + M 2 ¡ 1 y l 1 +1 ;l 2 y l 1 +1 ;l 2 +1 ¢ ¢ ¢ y l 1 +1 ;l 2 + M 2 ¡ 1 . . . . . . . . . . . . y l 1 + M 1 ¡ 1 ;l 2 y l 1 + M 1 ¡ 1 ;l 2 +1 ¢ ¢ ¢ y l 1 + M 1 ¡ 1 ;l 2 + M 2 ¡ 1 3 7 7 7 5 ; (4.16) where l 1 = 0 ; :::; L 1 ¡ 1 and l 2 = 0 ; :::; L 2 ¡ 1. Increasing M 1 and M 2 t ypically increases the sp ectral resolution at the cost of reducing the statistical stabilit y of the sp ectral estimates due to the reduced n um b er of submatrices. T ypically , w e c ho ose M 1 = N 1 = 2 and M 2 = N 2 = 2 [36 , 34]. Recall ¹ y l 1 ;l 2 = v ec[ ¹ Y l 1 ;l 2 ] (4.17) and a ( ! 1 ; ! 2 ) = a M 2 ( ! 2 ) ­ a M 1 ( ! 1 ) ; (4.18) where ­ denotes the Kronec k er matrix pro duct, and a M k ( ! k ) 4 = [1 e j ! k ¢ ¢ ¢ e j ( M k ¡ 1) ! k ] T ; k = 1 ; 2 : (4.19) Then, according to (4.1), the snapshot v ector ¹ y l 1 ;l 2 corresp onding to ¹ Y l 1 ;l 2 can b e written as ¹ y l 1 ;l 2 = [ ® ( ! 1 ; ! 2 ) a ( ! 1 ; ! 2 )] ¢ e j ( ! 1 l 1 + ! 2 l 2 ) + ¹ e l 1 ;l 2 ( ! 1 ; ! 2 ) ; (4.20) where ¹ e l 1 ;l 2 ( ! 1 ; ! 2 ) is formed from f e n 1 ;n 2 ( ! 1 ; ! 2 ) g in the same w a y as ¹ y l 1 ;l 2 is made from f y n 1 ;n 2 g . T o estimate ® ( ! 1 ; ! 2 ), the APES algorithm mimics an ML estimator b y assuming that f ¹ e l 1 ;l 2 ( ! 1 ; ! 2 ) g L 1 ¡ 1 ;L 2 ¡ 1 l 1 =0 ;l 2 =0 are zero-mean circularly symmetric complex Gaussian random v ectors that are statistically indep enden t of eac h other and ha v e

PAGE 64

52 the same unkno wn co v ariance matrix Q ( ! 1 ; ! 2 ) = E £ ¹ e l 1 ;l 2 ( ! 1 ; ! 2 ) ¹ e H l 1 ;l 2 ( ! 1 ; ! 2 ) ¤ : (4.21) Using the ab o v e assumptions, w e get the normalized surrogate log-lik eliho o d function of the data snapshots f ¹ y l 1 ;l 2 g as follo ws: 1 L 1 L 2 ln p ( f ¹ y l l ;l 2 g j ® ( ! 1 ; ! 2 ) ; Q ( ! 1 ; ! 2 ) ) = ¡ M 1 M 2 ln ¼ ¡ ln j Q ( ! 1 ; ! 2 ) j ¡ 1 L 1 L 2 L 1 ¡ 1 X l 1 =0 L 2 ¡ 1 X l 2 =0 [ ¹ y l 1 ;l 2 ¡ ® ( ! 1 ; ! 2 ) a ( ! 1 ; ! 2 ) e j ( ! 1 l 1 + ! 2 l 2 ) ¤ H Q ¡ 1 ( ! 1 ; ! 2 ) £ ¹ y l 1 ;l 2 ¡ ® ( ! 1 ; ! 2 ) a ( ! 1 ; ! 2 ) e j ( ! 1 l 1 + ! 2 l 2 ) ¤ (4.22) = ¡ M 1 M 2 ln ¼ ¡ ln j Q ( ! 1 ; ! 2 ) j ¡ tr ½ Q ¡ 1 ( ! 1 ; ! 2 ) ¢ 1 L 1 L 2 L 1 ¡ 1 X l 1 =0 L 2 ¡ 1 X l 2 =0 [ ¹ y l 1 ;l 2 ¡ ® ( ! 1 ; ! 2 ) a ( ! 1 ; ! 2 ) e j ( ! 1 l 1 + ! 2 l 2 ) ¤ £ ¹ y l 1 ;l 2 ¡ ® ( ! 1 ; ! 2 ) a ( ! 1 ; ! 2 ) e j ( ! 1 l 1 + ! 2 l 2 ) ¤ H o : (4.23) Similar to the 1-D case, the maximization of the ab o v e surrogate lik eliho o d function giv es the APES estimator ^ ® ( ! 1 ; ! 2 ) = a H ( ! 1 ; ! 2 ) ^ S ¡ 1 ( ! 1 ; ! 2 ) ¹ g ( ! 1 ; ! 2 ) a H ( ! 1 ; ! 2 ) ^ S ¡ 1 ( ! 1 ; ! 2 ) a ( ! 1 ; ! 2 ) (4.24) and ^ Q ( ! 1 ; ! 2 ) = ^ S ( ! 1 ; ! 2 ) + [ ^ ® M L ( ! 1 ; ! 2 ) a ( ! 1 ; ! 2 ) ¡ ¹ g ( ! 1 ; ! 2 )] [ ^ ® M L ( ! 1 ; ! 2 ) a ( ! 1 ; ! 2 ) ¡ ¹ g ( ! 1 ; ! 2 )] H ; (4.25) where ^ R , ¹ g ( ! 1 ; ! 2 ), and ^ S ( ! 1 ; ! 2 ) ha v e b een de¯ned in Section 4.2.

PAGE 65

CHAPTER 5 TW O-DIMENSIONAL MAPES VIA EXPECT A TION MAXIMIZA TION AND CYCLIC MAXIMIZA TION 5.1 In tro duction In Chapter 3, w e ha v e prop osed the one-dimensional (1-D) missing-data amplitude and phase estimation via exp ectation maximization (MAPES-EM) algorithms to deal with the general missing-data problem where the missing data samples o ccur in arbitrary patterns [80 ]. The MAPES-EM algorithms w ere deriv ed follo wing a \maxim um lik eliho o d (ML) ¯tting " based approac h in whic h a \ML ¯tting" problem w as solv ed iterativ ely via the exp ectation maximization (EM) algorithm. MAPES-EM exhibits b etter sp ectral estimation p erformance than GAPES. Ho w ev er, the MAPESEM algorithms are computationally in tensiv e due to estimating the missing samples separately for eac h frequency of in terest. The direct application of MAPES-EM to large data sets, e.g., t w o-dimensional (2-D) data, is computationally prohibitiv e. In this c hapter, w e consider the problem of 2-D nonparametric sp ectral estimation of data matrices with missing data samples o ccurring in arbitrary patterns [79 ]. First, w e presen t 2-D extensions of the MAPES-EM algorithms in tro duced in Chapter 3 in the 1-D case. Then w e dev elop a new MAPES algorithm, referred to as MAPES-CM, b y solving an ML problem iterativ ely via cyclic maximization (CM). MAPES-EM and MAPES-CM p ossess similar sp ectral estimation p erformance, but the computational complexit y of the latter is m uc h lo w er than that of the former. The remainder of this Chapter is organized as follo ws. In Section 5.2, w e presen t 2-D extensions of the MAPES-EM algorithms and dev elop the 2-D MAPESCM algorithm in Section 5.3. In Section 5.4, w e compare MAPES-CM with MAPES53

PAGE 66

54 EM, from b oth theoretical and computational p oin t-of-view. Numerical examples are pro vided in Section 5.5 to demonstrate the p erformance of the MAPES algorithms. Section 5.6 concludes this c hapter. 5.2 2-D MAPES via Exp ectation Maximization Assume that some arbitrary elemen ts of the data matrix Y are missing. Due to these missing data samples, whic h can b e treated as unkno wns, the log-lik eliho o d function (2.22) cannot b e maximized directly . In this section, w e will sho w ho w to tac kle this missing-data problem, in the ML con text, through the use of the exp ectation maximization (EM) and cyclic maximization (CM) algorithms. A comparison of these t w o approac hes is also pro vided. 5.2.1 2-D MAPES-EM1 W e assume that the data snapshots f ¹ Y l 1 ;l 2 g (or f ¹ y l 1 ;l 2 g ) are indep enden t of eac h other, and w e estimate the missing data separately for di®eren t data snapshots. F or eac h data snapshot ¹ y l 1 ;l 2 , let ¹ ° l 1 ;l 2 and ¹ ¹ l 1 ;l 2 denote the v ectors con taining the a v ailable and missing elemen ts of ¹ y l 1 ;l 2 , resp ectiv ely . Assume that ¹ ° l 1 ;l 2 has dimension g l 1 ;l 2 £ 1 where 1 · g l 1 ;l 2 · M 1 M 2 is the n um b er of a v ailable elemen ts in the snapshot ¹ y l 1 ;l 2 . Then ¹ ° l 1 ;l 2 and ¹ ¹ l 1 ;l 2 are related to ¹ y l 1 ;l 2 b y unitary transformations as follo ws: ¹ ° l 1 ;l 2 = ¹ S T g ( l 1 ; l 2 ) ¹ y l 1 ;l 2 (5.1) ¹ ¹ l 1 ;l 2 = ¹ S T m ( l 1 ; l 2 ) ¹ y l 1 ;l 2 ; (5.2) where ¹ S g ( l 1 ; l 2 ) and ¹ S m ( l 1 ; l 2 ) are M 1 M 2 £ g l 1 ;l 2 and M 1 M 2 £ ( M 1 M 2 ¡ g l 1 ;l 2 ) unitary selection matrices suc h that ¹ S T g ( l 1 ; l 2 ) ¹ S g ( l 1 ; l 2 ) = I g l 1 ;l 2 , ¹ S T m ( l 1 ; l 2 ) ¹ S m ( l 1 ; l 2 ) = I M 1 M 2 ¡ g l 1 ;l 2 and ¹ S T g ( l 1 ; l 2 ) ¹ S m ( l 1 ; l 2 ) = 0 g l 1 ;l 2 £ ( M 1 M 2 ¡ g l 1 ;l 2 ) . F or example, let M 1 = 3, M 2 = 2, and let ¹ Y l 1 ;l 2 = 2 4 y l 1 ;l 2 ? ? ? y l 1 +2 ;l 2 y l 1 +2 ;l 2 +1 3 5 ; (5.3)

PAGE 67

55 where eac h ? indicates a missing sample. Then w e ha v e g l 1 ;l 2 = 3, ¹ y l 1 ;l 2 = [ y l 1 ;l 2 ? y l 1 +2 ;l 2 ? ? y l 1 +2 ;l 2 +1 ] T ; (5.4) and ¹ S g ( l 1 ; l 2 ) = 2 6 6 6 6 6 6 4 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 3 7 7 7 7 7 7 5 ; ¹ S m ( l 1 ; l 2 ) = 2 6 6 6 6 6 6 4 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 0 0 0 3 7 7 7 7 7 7 5 : (5.5) Because w e ha v e ¹ y l 1 ;l 2 = £ ¹ S g ( l 1 ; l 2 ) ¹ S T g ( l 1 ; l 2 ) + ¹ S m ( l 1 ; l 2 ) ¹ S T m ( l 1 ; l 2 ) ¤ ¹ y l 1 ;l 2 = ¹ S g ( l 1 ; l 2 ) ¹ ° l 1 ;l 2 + ¹ S m ( l 1 ; l 2 ) ¹ ¹ l 1 ;l 2 ; (5.6) the join t normalized surrogate log-lik eliho o d function of f ¹ ° l 1 ;l 2 ; ¹ ¹ l 1 ;l 2 g is obtained b y substituting (5.6) in to (4.23) 1 L 1 L 2 ln p ( f ¹ ° l 1 ;l 2 ; ¹ ¹ l 1 ;l 2 g j ® ( ! 1 ; ! 2 ) ; Q ( ! 1 ; ! 2 ) ) = ¡ M 1 M 2 ln ¼ ¡ ln j Q ( ! 1 ; ! 2 ) j ¡ tr ½ Q ¡ 1 ( ! 1 ; ! 2 ) ¢ 1 L 1 L 2 L 1 ¡ 1 X l 1 =0 L 2 ¡ 1 X l 2 =0 £ ¹ S g ( l 1 ; l 2 ) ¹ ° l 1 ;l 2 + ¹ S m ( l 1 ; l 2 ) ¹ ¹ l 1 ;l 2 ¡ ® ( ! 1 ; ! 2 ) a ( ! 1 ; ! 2 ) e j ( ! 1 l 1 + ! 2 l 2 ) ¤ £ ¹ S g ( l 1 ; l 2 ) ¹ ° l 1 ;l 2 + ¹ S m ( l 1 ; l 2 ) ¹ ¹ l 1 ;l 2 ¡ ® ( ! 1 ; ! 2 ) a ( ! 1 ; ! 2 ) e j ( ! 1 l 1 + ! 2 l 2 ) ¤ H o : (5.7) Similar to the 1-D case, the probabilit y densit y function of ¹ ¹ l 1 ;l 2 conditioned on ¹ ° l 1 ;l 2 (for giv en µ = ^ µ i ¡ 1 ) is complex Gaussian with mean ¹ b l 1 ;l 2 and co v ariance matrix ¹ K l 1 ;l 2 : ¹ ¹ l 1 ;l 2 j ¹ ° l 1 ;l 2 ; ^ µ i ¡ 1 » C N ( ¹ b l 1 ;l 2 ; ¹ K l 1 ;l 2 ) ; (5.8) where ¹ b l 1 ;l 2 = E n ¹ ¹ l 1 ;l 2 ¯ ¯ ¯ ¹ ° l 1 ;l 2 ; ^ µ i ¡ 1 o = ¹ S T m ( l 1 ; l 2 ) a ( ! 1 ; ! 2 ) ^ ® i ¡ 1 ( ! 1 ; ! 2 ) e j ( ! 1 l 1 + ! 2 l 2 ) + ¹ S T m ( l 1 ; l 2 ) ^ Q i ¡ 1 ( ! 1 ; ! 2 ) ¹ S g ( l 1 ; l 2 ) h ¹ S T g ( l 1 ; l 2 ) ^ Q i ¡ 1 ( ! 1 ; ! 2 ) ¹ S g ( l 1 ; l 2 ) i ¡ 1 £ ¹ ° l 1 ;l 2 ¡ ¹ S T g ( l 1 ; l 2 ) a ( ! 1 ; ! 2 ) ^ ® i ¡ 1 ( ! 1 ; ! 2 ) e j ( ! 1 l 1 + ! 2 l 2 ) ¤ (5.9)

PAGE 68

56 and ¹ K l 1 ;l 2 = co v n ¹ ¹ l 1 ;l 2 ¯ ¯ ¯ ¹ ° l 1 ;l 2 ; ^ µ i ¡ 1 o = ¹ S T m ( l 1 ; l 2 ) ^ Q i ¡ 1 ( ! 1 ; ! 2 ) ¹ S m ( l 1 ; l 2 ) ¡ ¹ S T m ( l 1 ; l 2 ) ^ Q i ¡ 1 ( ! 1 ; ! 2 ) ¹ S g ( l 1 ; l 2 ) h ¹ S T g ( l 1 ; l 2 ) ^ Q i ¡ 1 ( ! 1 ; ! 2 ) ¹ S g ( l 1 ; l 2 ) i ¡ 1 ¹ S T g ( l 1 ; l 2 ) ^ Q i ¡ 1 ( ! 1 ; ! 2 ) ¹ S m ( l 1 ; l 2 ) : (5.10) 5.2.1.1 Exp ectation The conditional exp ectation of the surrogate log-lik eliho o d in (5.7) according to (5.8)-(5.10) is giv en b y E ½ 1 L 1 L 2 ln p ( © ¹ ° l 1 ;l 2 ; ¹ ¹ l 1 ;l 2 ª j ® ( ! 1 ; ! 2 ) ; Q ( ! 1 ; ! 2 ) ) ¯ ¯ © ¹ ° l 1 ;l 2 ª ; ^ ® i ¡ 1 ( ! 1 ; ! 2 ) ; ^ Q i ¡ 1 ( ! 1 ; ! 2 ) o = ¡ M 1 M 2 ln ¼ ¡ ln j Q ( ! 1 ; ! 2 ) j ¡ tr ( Q ¡ 1 ( ! 1 ; ! 2 ) 1 L 1 L 2 L 1 ¡ 1 X l 1 =0 L 2 ¡ 1 X l 2 =0 Ã ¹ S m ( l 1 ; l 2 ) ¹ K l 1 ;l 2 ¹ S T m ( l 1 ; l 2 ) + £ ¹ S g ( l 1 ; l 2 ) ¹ ° l 1 ;l 2 + ¹ S m ( l 1 ; l 2 ) ¹ b l 1 ;l 2 ¡ ® ( ! 1 ; ! 2 ) a ( ! 1 ; ! 2 ) e j ( ! 1 l 1 + ! 2 l 2 ) ¤ £ ¹ S g ( l 1 ; l 2 ) ¹ ° l 1 ;l 2 + ¹ S m ( l 1 ; l 2 ) ¹ b l 1 ;l 2 ¡ ® ( ! 1 ; ! 2 ) a ( ! 1 ; ! 2 ) e j ( ! 1 l 1 + ! 2 l 2 ) ¤ H ´ o : (5.11) 5.2.1.2 Maximization The normalized exp ected surrogate log-lik eliho o d (5.11) can b e re-written as ¡ M 1 M 2 ln ¼ ¡ ln j Q ( ! 1 ; ! 2 ) j ¡ tr ( Q ¡ 1 ( ! 1 ; ! 2 ) 1 L 1 L 2 L 1 ¡ 1 X l 1 =0 L 2 ¡ 1 X l 2 =0 Ã ¹ ¡ l 1 ;l 2 + [ ¹ z l 1 ;l 2 ¡ ® ( ! 1 ; ! 2 ) a ( ! 1 ; ! 2 ) e j ( ! 1 l 1 + ! 2 l 2 ) ¤ £ ¹ z l 1 ;l 2 ¡ ® ( ! 1 ; ! 2 ) a ( ! 1 ; ! 2 ) e j ( ! 1 l 1 + ! 2 l 2 ) ¤ H ´ o ; (5.12) where w e ha v e de¯ned ¹ ¡ l 1 ;l 2 4 = ¹ S m ( l 1 ; l 2 ) ¹ K l 1 ;l 2 ¹ S T m ( l 1 ; l 2 ) (5.13) and ¹ z l 1 ;l 2 4 = ¹ S g ( l 1 ; l 2 ) ¹ ° l 1 ;l 2 + ¹ S m ( l 1 ; l 2 ) ¹ b l 1 ;l 2 : (5.14) Maximizing (5.12) giv es ^ ® 1 ( ! 1 ; ! 2 ) = a H ( ! 1 ; ! 2 ) ¹ S ¡ 1 ( ! 1 ; ! 2 ) ¹ Z ( ! 1 ; ! 2 ) a H ( ! 1 ; ! 2 ) ¹ S ¡ 1 ( ! 1 ; ! 2 ) a ( ! 1 ; ! 2 ) (5.15)

PAGE 69

57 and ^ Q 1 ( ! 1 ; ! 2 ) = ¹ S ( ! 1 ; ! 2 ) + £ ^ ® 1 ( ! 1 ; ! 2 ) a ( ! 1 ; ! 2 ) ¡ ¹ Z ( ! 1 ; ! 2 ) ¤ £ ^ ® 1 ( ! 1 ; ! 2 ) a ( ! 1 ; ! 2 ) ¡ ¹ Z ( ! 1 ; ! 2 ) ¤ H ; (5.16) where ¹ Z ( ! 1 ; ! 2 ) 4 = 1 L 1 L 2 L 1 ¡ 1 X l 1 =0 L 2 ¡ 1 X l 2 =0 ¹ z l 1 ;l 2 e ¡ j ( ! 1 l 1 + ! 2 l 2 ) (5.17) and ¹ S ( ! 1 ; ! 2 ) 4 = 1 L 1 L 2 L 1 ¡ 1 X l 1 =0 L 2 ¡ 1 X l 2 =0 ¹ ¡ l 1 ;l 2 + 1 L 1 L 2 L 1 ¡ 1 X l 1 =0 L 2 ¡ 1 X l 2 =0 ¹ z l 1 ;l 2 ¹ z H l 1 ;l 2 ¡ ¹ Z ( ! 1 ; ! 2 ) ¹ Z H ( ! 1 ; ! 2 ) : (5.18) This completes the deriv ation of the 2-D MAPES-EM1 algorithm, a step-b y-step summary of whic h is giv en b elo w. Step 0 : Obtain an initial estimate of f ® ( ! 1 ; ! 2 ) ; Q ( ! 1 ; ! 2 ) g . Step 1 : Use the most recen t estimate of f ® ( ! 1 ; ! 2 ) ; Q ( ! 1 ; ! 2 ) g in (5.9) and (5.10) to calculate ¹ b l 1 ;l 2 and ¹ K l 1 ;l 2 , resp ectiv ely . Note that ¹ b l 1 ;l 2 can b e regarded as the curren t estimate of the corresp onding missing samples. Step 2 : Up date the estimate of f ® ( ! 1 ; ! 2 ) ; Q ( ! 1 ; ! 2 ) g using (5.15) and (5.16). Step 3 : Rep eat Steps 1 and 2 un til practical con v ergence. 5.2.2 2-D MAPES-EM2 MAPES-EM2 utilizes the EM algorithm b y estimating the missing data sim ultaneously for all data snapshots. Let y = v ec [ Y ] (5.19) denotes the v ector of all the data samples. Recall that ° and ¹ denote the v ectors con taining the a v ailable and missing elemen ts of y , resp ectiv ely , where ° has a size of g £ 1.

PAGE 70

58 W e let ¸ y denote the L 1 L 2 M 1 M 2 £ 1 v ector obtained b y concatenating all the snapshots ¸ y 4 = 2 6 4 ¹ y 0 ; 0 . . . ¹ y L 1 ¡ 1 ;L 2 ¡ 1 3 7 5 = S g ° + S m ¹ ; (5.20) where S g (whic h has a size of L 1 L 2 M 1 M 2 £ g ) and S m (whic h has a size of L 1 L 2 M 1 M 2 £ ( N 1 N 2 ¡ g )) are the corresp onding selection matrices for the a v ailable and missing data v ectors, resp ectiv ely . Due to the o v erlapping of the v ectors f ¹ y l 1 ;l 2 g , S g and S m are not unitary , but they are still orthogonal to eac h other: S T g S m = 0 g £ ( N 1 N 2 ¡ g ) . So instead of (5.1) and (5.2), w e ha v e from (5.20): ° = ( S T g S g ) ¡ 1 S T g ¸ y = ¸ S T g ¸ y (5.21) and ¹ = ( S T m S m ) ¡ 1 S T m ¸ y = ¸ S T m ¸ y ; (5.22) where the matrices ¸ S g and ¸ S m in tro duced ab o v e are de¯ned as ¸ S g 4 = S g ( S T g S g ) ¡ 1 , ¸ S m 4 = S m ( S T m S m ) ¡ 1 , and they are also orthogonal to eac h other: ¸ S T g ¸ S m = 0 g £ ( N 1 N 2 ¡ g ) . No w the normalized surrogate log-lik eliho o d function in (4.22) can b e written as 1 L 1 L 2 ln p ( ¸ y j ® ( ! 1 ; ! 2 ) ; Q ( ! 1 ; ! 2 ) ) = ¡ M 1 M 2 ln ¼ ¡ 1 L 1 L 2 ln j D ( ! 1 ; ! 2 ) j ¡ 1 L 1 L 2 [ ¸ y ¡ ® ( ! 1 ; ! 2 ) ½ ( ! 1 ; ! 2 )] H D ¡ 1 ( ! 1 ; ! 2 ) [ ¸ y ¡ ® ( ! 1 ; ! 2 ) ½ ( ! 1 ; ! 2 )] ; (5.23) where ½ ( ! 1 ; ! 2 ) and D ( ! 1 ; ! 2 ) are de¯ned as: ½ ( ! 1 ; ! 2 ) 4 = 2 6 4 e j ( ! 1 0+ ! 2 0) a ( ! 1 ; ! 2 ) . . . e j [ ! 1 ( L 1 ¡ 1)+ ! 2 ( L 2 ¡ 1)] a ( ! 1 ; ! 2 ) 3 7 5 (5.24) and D ( ! 1 ; ! 2 ) 4 = 2 6 4 Q ( ! 1 ; ! 2 ) 0 . . . 0 Q ( ! 1 ; ! 2 ) 3 7 5 : (5.25)

PAGE 71

59 Substituting (5.20) in to (5.23), w e obtain the join t surrogate log-lik eliho o d of ° and ¹ : 1 L 1 L 2 ln p ( ° ; ¹ j ® ( ! 1 ; ! 2 ) ; Q ( ! 1 ; ! 2 ) ) = 1 L 1 L 2 ½ ¡ L 1 L 2 M 1 M 2 ln ¼ ¡ ln j D ( ! 1 ; ! 2 ) j ¡ [ S g ° + S m ¹ ¡ ® ( ! 1 ; ! 2 ) ½ ( ! 1 ; ! 2 )] H D ¡ 1 ( ! 1 ; ! 2 ) [ S g ° + S m ¹ ¡ ® ( ! 1 ; ! 2 ) ½ ( ! 1 ; ! 2 )] o + C J ; (5.26) where C J is a constan t that accoun ts for the Jacobian of the non unitary transformation b et w een ¸ y and ° and ¹ in (5.20). T o deriv e the EM algorithm for the curren t set of assumptions, w e ha v e ¹ j ° ; ^ µ i ¡ 1 » C N ( b ; K ) (5.27) where b = E n ¹ ¯ ¯ ¯ ° ; ^ µ i ¡ 1 o = ¸ S T m ½ ( ! 1 ; ! 2 ) ^ ® i ¡ 1 ( ! 1 ; ! 2 ) + ¸ S T m ^ D i ¡ 1 ( ! 1 ; ! 2 ) ¸ S g h ¸ S T g ^ D i ¡ 1 ( ! 1 ; ! 2 ) ¸ S g i ¡ 1 h ° ¡ ¸ S T g ½ ( ! 1 ; ! 2 ) ^ ® i ¡ 1 ( ! 1 ; ! 2 ) i (5.28) and K = co v n ¹ ¯ ¯ ¯ ° ; ^ µ i ¡ 1 o = ¸ S T m ^ D i ¡ 1 ( ! 1 ; ! 2 ) ¸ S m ¡ ¸ S T m ^ D i ¡ 1 ( ! 1 ; ! 2 ) ¸ S g h ¸ S T g ^ D i ¡ 1 ( ! 1 ; ! 2 ) ¸ S g i ¡ 1 ¸ S T g ^ D i ¡ 1 ( ! 1 ; ! 2 ) ¸ S m : (5.29) 5.2.2.1 Exp ectation The conditional exp ectation of the surrogate log-lik eliho o d function in (5.26) is giv en as E ½ 1 L 1 L 2 ln p ( ° ; ¹ j ® ( ! 1 ; ! 2 ) ; Q ( ! 1 ; ! 2 ) ) ¯ ¯ ¯ ° ; ^ ® i ¡ 1 ( ! 1 ; ! 2 ) ; ^ Q i ¡ 1 ( ! 1 ; ! 2 ) ¾ = ¡ M 1 M 2 ln ¼ ¡ 1 L 1 L 2 ln j D ( ! 1 ; ! 2 ) j ¡ tr ½ 1 L 1 L 2 D ¡ 1 ( ! 1 ; ! 2 ) µ S m KS T m + [ S g ° + S m b ¡ ® ( ! 1 ; ! 2 ) ½ ( ! 1 ; ! 2 )] [ S g ° + S m b ¡ ® ( ! 1 ; ! 2 ) ½ ( ! 1 ; ! 2 )] H ´o + C J : (5.30)

PAGE 72

60 5.2.2.2 Maximization T o maximize the exp ected surrogate log-lik eliho o d function in (5.30), w e again exploit the kno wn structure of D ( ! 1 ; ! 2 ) and ½ ( ! 1 ; ! 2 ). Let 2 6 4 z 0 ; 0 . . . z L 1 ¡ 1 ;L 2 ¡ 1 3 7 5 4 = S g ° + S m b (5.31) denote the data snapshots made up of the a v ailable and estimated data samples, where eac h z l 1 ;l 2 ( l 1 = 0 ; :::; L 1 ¡ 1 and l 2 = 0 ; :::; L 2 ¡ 1) is an M 1 M 2 £ 1 v ector. Also let ¡ 0 ; 0 , ..., ¡ L 1 ¡ 1 ;L 2 ¡ 1 b e the M 1 M 2 £ M 1 M 2 blo c ks on the blo c k diagonal of S m KS T m . Then the exp ected surrogate log-lik eliho o d function w e need to maximize with resp ect to ® ( ! 1 ; ! 2 ) and Q ( ! 1 ; ! 2 ) b ecomes (to within an additiv e constan t) ¡ ln j Q ( ! 1 ; ! 2 ) j ¡ tr ½ Q ¡ 1 ( ! 1 ; ! 2 ) 1 L 1 L 2 L 1 ¡ 1 X l 1 =0 L 2 ¡ 1 X l 2 =0 ( ¡ l 1 ;l 2 + [ z l 1 ;l 2 ¡ ® ( ! 1 ; ! 2 ) a ( ! 1 ; ! 2 ) e j ( ! 1 l 1 + ! 2 l 2 ) ¤ £ z l 1 ;l 2 ¡ ® ( ! 1 ; ! 2 ) a ( ! 1 ; ! 2 ) e j ( ! 1 l 1 + ! 2 l 2 ) ¤ H ´o : (5.32) The solution b ecomes ^ ® 2 ( ! 1 ; ! 2 ) = a H ( ! 1 ; ! 2 ) S ¡ 1 ( ! 1 ; ! 2 ) Z ( ! 1 ; ! 2 ) a H ( ! 1 ; ! 2 ) S ¡ 1 ( ! 1 ; ! 2 ) a ( ! 1 ; ! 2 ) (5.33) and ^ Q 2 ( ! 1 ; ! 2 ) = S ( ! 1 ; ! 2 ) + [ ^ ® 2 ( ! 1 ; ! 2 ) a ( ! 1 ; ! 2 ) ¡ Z ( ! 1 ; ! 2 )] [ ^ ® 2 ( ! 1 ; ! 2 ) a ( ! 1 ; ! 2 ) ¡ Z ( ! 1 ; ! 2 )] H ; (5.34) where S ( ! 1 ; ! 2 ) and Z ( ! 1 ; ! 2 ) are de¯ned as S ( ! 1 ; ! 2 ) 4 = 1 L 1 L 2 L 1 ¡ 1 X l 1 =0 L 2 ¡ 1 X l 2 =0 ¡ l 1 ;l 2 + 1 L 1 L 2 L 1 ¡ 1 X l 1 =0 L 2 ¡ 1 X l 2 =0 z l 1 ;l 2 z H l 1 ;l 2 ¡ Z ( ! 1 ; ! 2 ) Z H ( ! 1 ; ! 2 ) (5.35) and Z ( ! 1 ; ! 2 ) 4 = 1 L 1 L 2 L 1 ¡ 1 X l 1 =0 L 2 ¡ 1 X l 2 =0 z l 1 ;l 2 e ¡ j ( ! 1 l 1 + ! 2 l 2 ) : (5.36) The deriv ation of the MAPES-EM2 algorithm is th us complete, and a step-b y-step summary of this algorithm follo ws:

PAGE 73

61 Step 0 : Obtain an initial estimate of f ® ( ! 1 ; ! 2 ) ; Q ( ! 1 ; ! 2 ) g . Step 1 : Use the most recen t estimates of f ® ( ! 1 ; ! 2 ) ; Q ( ! 1 ; ! 2 ) g in (5.28) and (5.29) to calculate b and K . Note that b can b e regarded as the curren t estimate of the missing sample v ector. Step 2 : Up date the estimates of f ® ( ! 1 ; ! 2 ) ; Q ( ! 1 ; ! 2 ) g using (5.33) and (5.34). Step 3 : Rep eat Steps 1 and 2 un til practical con v ergence. 5.3 2-D MAPES via Cyclic Maximization Next w e consider ev aluating the sp ectrum on the K 1 £ K 2 -p oin t DFT grid. Instead of dealing with eac h individual frequency ( ! k 1 ; ! k 2 ) separately , w e consider the follo wing maximization problem: max ¹ ; f ® ( ! k 1 ;! k 2 ) ; Q ( ! k 1 ;! k 2 ) g K 1 ¡ 1 X k 1 =0 K 2 ¡ 1 X k 2 =0 ( ¡ ln j Q ( ! k 1 ; ! k 2 ) j ¡ 1 L 1 L 2 L 1 ¡ 1 X l 1 =0 L 2 ¡ 1 X l 2 =0 £ ¹ y l 1 ;l 2 ¡ ® ( ! k 1 ; ! k 2 ) a ( ! k 1 ; ! k 2 ) e j ( ! k 1 l 1 + ! k 2 l 2 ) ¤ H Q ¡ 1 ( ! k 1 ; ! k 2 ) £ ¹ y l 1 ;l 2 ¡ ® ( ! k 1 ; ! k 2 ) a ( ! k 1 ; ! k 2 ) e j ( ! k 1 l 1 + ! k 2 l 2 ) ¤ ª ; (5.37) where the ob jectiv e function is the summation o v er the 2-D frequency grid of all the frequency-dep enden t complete-data lik eliho o d functions in (4.22) (within an additiv e constan t). W e solv e the ab o v e optimization problem via a cyclic maximization (CM) approac h. First, assuming that the previous estimate ^ µ i ¡ 1 formed from f ^ ® i ¡ 1 ( ! k 1 ; ! k 2 ) ; ^ Q i ¡ 1 ( ! k 1 ; ! k 2 ) g is a v ailable, w e maximize (5.37) with resp ect to ¹ . This step can b e re-form ulated as min ¹ K 1 ¡ 1 X k 1 =0 K 2 ¡ 1 X k 2 =0 £ ¸ y ¡ ^ ® i ¡ 1 ( ! k 1 ; ! k 2 ) ½ ( ! k 1 ; ! k 2 ) ¤ H h ^ D i ¡ 1 ( ! k 1 ; ! k 2 ) i ¡ 1 [ ¸ y ¡ ^ ® i ¡ 1 ( ! k 1 ; ! k 2 ) ½ ( ! k 1 ; ! k 2 ) ¤ ; (5.38) where ¸ y , ½ ( ! k 1 ; ! k 2 ), and ^ D i ¡ 1 ( ! 1 ; ! 2 ) ha v e b een de¯ned previously . Recalling that ¸ y = S g ° + S m ¹ , w e can easily solv e the optimization problem in (5.38) as its ob jectiv e

PAGE 74

62 function is quadratic in ¹ : ^ ¹ = " K 1 ¡ 1 X k 1 =0 K 2 ¡ 1 X k 2 =0 S H m h ^ D i ¡ 1 ( ! k 1 ; ! k 2 ) i ¡ 1 S m # ¡ 1 K 1 ¡ 1 X k 1 =0 K 2 ¡ 1 X k 2 =0 S H m h ^ D i ¡ 1 ( ! k 1 ; ! k 2 ) i ¡ 1 £ ^ ® i ¡ 1 ( ! k 1 ; ! k 2 ) ½ ( ! k 1 ; ! k 2 ) ¡ S g ° ¤ : (5.39) A necessary condition for the in v erse in (5.39) to exist is that L 1 L 2 M 1 M 2 > N 1 N 2 ¡ g , whic h is alw a ys satis¯ed. Once an estimate ^ ¹ has b ecome a v ailable, w e re-estimate f ® ( ! k 1 ; ! k 2 ) g and f Q ( ! k 1 ; ! k 2 ) g b y maximizing (5.37) with ¹ replaced b y ^ ¹ . This can b e done b y maximizing eac h frequency term separately: max ® ( ! k 1 ;! k 2 ) ; Q ( ! k 1 ;! k 2 ) ¡ ln j Q ( ! k 1 ; ! k 2 ) j ¡ 1 L 1 L 2 L 1 ¡ 1 X l 1 =0 L 2 ¡ 1 X l 2 =0 £ ^ ¹ y l 1 ;l 2 ¡ ® ( ! k 1 ; ! k 2 ) a ( ! k 1 ; ! k 2 ) e j ( ! k 1 l 1 + ! k 2 l 2 ) ¤ H Q ¡ 1 ( ! k 1 ; ! k 2 ) £ ^ ¹ y l 1 ;l 2 ¡ ® ( ! k 1 ; ! k 2 ) a ( ! k 1 ; ! k 2 ) e j ( ! k 1 l 1 + ! k 2 l 2 ) ¤ ; (5.40) whic h reduces to the 2-D APES problem. A cyclic maximization of (5.37) can b e implemen ted b y the alternating maximization with resp ect to ¹ and, resp ectiv ely , ® ( ! k 1 ; ! k 2 ) and Q ( ! k 1 ; ! k 2 ). A step-b ystep summary of 2-D MAPES-CM is as follo ws. Step 0 : Obtain an initial estimate of f ® ( ! k 1 ; ! k 2 ) ; Q ( ! k 1 ; ! k 2 ) g . Step 1 : Use the most recen t estimates of f ® ( ! k 1 ; ! k 2 ) ; Q ( ! k 1 ; ! k 2 ) g in (5.39) to estimate the missing samples. Step 2 : Up date the estimates of f ® ( ! 1 ; ! 2 ) ; Q ( ! k 1 ; ! k 2 ) g using 2-D APES applied to the data matrices with the missing sample estimates from Step 1 (see (4.24) and (4.25)). Step 3 : Rep eat Steps 1 and 2 un til practical con v ergence.

PAGE 75

63 5.4 MAPES-EM v ersus MAPES-CM 5.4.1 P erformance Analysis Starting from the same ML based approac h, MAPES-EM and MAPES-CM end up di®eren tly b y maximizing their o wn lik eliho o d based ob jectiv e functions via di®eren t approac hes: MAPES-CM deals with all the frequencies sim ultaneously , and tries to maximize the summation (see (5.37)) of the frequency-dep enden t loglik eliho o d functions o v er a ¯xed DFT grid; MAPES-EM deals with eac h frequency separately , and tries to maximize the conditional exp ectation (see (5.11) and (5.30)) of the complete-data log-lik eliho o d function at eac h iteration. Since b oth MAPES-EM and MAPES-CM can b e expressed as a recursion ^ µ i = F ( ^ µ i ¡ 1 ) ; (5.41) a theoretical comparison of these t w o approac hes is made p ossible b y in v estigating their di®eren t recursiv e functions F , whic h corresp ond to Steps 1 and 2 in our stepb y-step summaries. ² Based on ^ µ i ¡ 1 , Step 1 of the MAPES algorithms tries to estimate the missing samples (e.g., ^ ¹ , as in MAPES-CM) or their statistics (e.g., ¹ b l 1 ;l 2 , b , ^ K l 1 ;l 2 , and K , as in MAPES-EM). The fundamen tal di®erence b et w een these t w o metho ds is that MAPES-EM utilizes the kno wn (a v ailable) data and further up dates its data mo del b y taking a Ba y esian approac h. This can b e understo o d b y considering a simple example. Assume that the distribution of a random v ector x is kno wn: x » C N ( b x ; K x ) : (5.42) If no comp onen t of x has b een observ ed, then w e can estimate x as ^ x M L = b x ; (5.43) whic h is also a minim um mean-square-error (MMSE) estimator. If some comp onen ts of x are a v ailable, sa y ° x , then t w o di®eren t approac hes can b e used to

PAGE 76

64 estimate the unkno wn comp onen ts ¹ x . One is to substitute the a v ailable data ° x in to the original ML problem and get the estimate of ¹ x b y maximizing the so-obtained lik eliho o d function: ^ ¹ x = arg min ¹ x ( x ¡ b x ) H K ¡ 1 x ( x ¡ b x ) : (5.44) The other is to consider the distribution of ¹ x conditioned on ° x : ¹ x j ° x » C N ( b ¹ j ° ; K ¹ j ° ) ; (5.45) and estimate ¹ x as ^ ¹ x = b ¹ j ° : (5.46) Clearly , (5.46) is a b etter estimate compared with (5.44), due to the fact that the former giv es the MMSE for the giv en ° x . By incorp orating a Ba y esian approac h, w e exp ect that MAPES-EM will pro vide more accurate sp ectral estimates than MAPES-CM. ² Step 2 of MAPES is to up date the estimates of f ® ( ! 1 ; ! 2 ) ; Q ( ! 1 ; ! 2 ) g according to (5.15)-(5.16) and (5.33)-(5.34). Comparing (5.15) and (5.33) with (4.24), w e can see that MAPES estimates the missing data and then uses these estimates ¹ b l 1 ;l 2 , b , or ^ ¹ as though they w ere correct. Ho w ev er, for the unstructured estimate ^ S ( ! 1 ; ! 2 ), MAPES-EM adds an extra term (see (5.18) or (5.35)) in v olving the conditional co v ariance ¹ K l 1 ;l 2 (or K ), whic h can b e regarded as a generalized diagonal loading op eration to mak e the sp ectral estimate robust against estimation errors. This loading op eration suggests that MAPES-EM will ha v e a slo w con v ergence. Comparing (5.16) and (5.34) with (4.25), the same loading op eration can b e observ ed when up dating the ML estimates of Q ( ! 1 ; ! 2 ). 5.4.2 Computational Analysis Consider ev aluating the sp ectrum for all three MAPES algorithms on the same DFT grid. Since all three algorithms iterate Step 1 and Step 2 un til practical

PAGE 77

65 con v ergence, w e can compare their computational complexit y separately for eac h step. ² In Step 1, MAPES-CM estimates the missing samples via (5.39), whic h can b e re-written in a simpli¯ed form as ^ ¹ = £ S H m D s S m ¤ ¡ 1 £ S H m D v ¡ S H m D s S g ° ¤ ; (5.47) where D s 4 = K 1 ¡ 1 X k 1 =0 K 2 ¡ 1 X k 2 =0 h ^ D i ¡ 1 ( ! k 1 ; ! k 2 ) i ¡ 1 (5.48) and D v 4 = K 1 ¡ 1 X k 1 =0 K 2 ¡ 1 X k 2 =0 h ^ D i ¡ 1 ( ! k 1 ; ! k 2 ) i ¡ 1 ^ ® i ¡ 1 ( ! k 1 ; ! k 2 ) ½ ( ! k 1 ; ! k 2 ) : (5.49) When computing D s and D v , the fact that ^ D i ¡ 1 ( ! k 1 ; ! k 2 ) is blo c k diagonal can b e exploited to reduce the computational complexit y . Comparing (5.47) with (5.28) and (5.29) (or (5.9) and (5.10)), whic h ha v e to b e ev aluated for eac h frequency ( ! k 1 ; ! k 2 ) (and for eac h snapshot ( l 1 ; l 2 ) for (5.9) and (5.10)), w e note that the computational complexit y of MAPES-CM is m uc h lo w er. ² In Step 2, MAPES-CM uses the standard APES algorithm whic h can b e e±cien tly implemen ted [39 , 31 ] as discussed in Section 2.6. As for the MAPES-EM sp ectral estimators ^ ® 1 ( ! 1 ; ! 2 ) and ^ ® 2 ( ! 1 ; ! 2 ), they ha v e the same structure as the APES estimator in (4.24), but with t w o di®erences. The ¯rst one is using the conditional mean ( ¹ b l 1 ;l 2 or b ) of the missing samples as their estimates. The second one is the additional terms ( ¹ ¡ l 1 ;l 2 and ¡ l 1 ;l 2 ) whic h in v olv e the cov ariance matrices ( ¹ K l 1 ;l 2 and K ) of the missing data in the ¹ S and S matrices. Due to the fact that MAPES-EM uses di®eren t estimates for the missing data at di®eren t frequencies, those tec hniques used to e±cien tly implemen t APES cannot b e applied here (see Section 2.7). As a result, no e±cien t algorithms are a v ailable to calculate the corresp onding sp ectral estimate.

PAGE 78

66 In summary , the MAPES-CM algorithm p ossesses a computational complexit y that is m uc h lo w er than that of MAPES-EM. 5.5 Numerical Examples In this section, w e presen t three n umerical examples to illustrate the p erformance of the MAPES algorithms for the 2-D missing-data sp ectral estimation problem. W e compare the MAPES algorithms with the windo w ed FFT (WFFT) and the GAPES [32 ]. A T a ylor windo w with order 5 and sidelob e lev el -35 dB is used to obtain the WFFT sp ectral estimates. All the 2-D sp ectra are plotted with a dynamic range of 35 dB. Similar to the 1-D case, w e simply let the initial estimate of ® ( ! 1 ; ! 2 ) b e giv en b y the windo w ed FFT (WFFT) with the missing data samples set to zero. The initial estimate of Q ( ! 1 ; ! 2 ) follo ws from the 2-D coun terpart of (2.24), where again, the missing data samples are set to zero. F or the initialization of GAPES, w e consider t w o cases. If the data missing pattern is arbitrary and no initial ¯lter with a prop er size can b e c hosen, w e use the same initial estimates of ® ( ! 1 ; ! 2 ) and Q ( ! 1 ; ! 2 ) as for MAPES, and the initial estimate of H ( ! 1 ; ! 2 ) follo ws from (4.13), where the missing samples are set to zero. If the data is gapp ed, as in Sections 5.5.2 and 5.5.3, w e follo w the initialization step in App endix 5.7. W e stop the iteration of all iterativ e algorithms whenev er the relativ e c hange in the total p o w er of the 2-D sp ectra corresp onding to the curren t ( ^ ® i ( ! k 1 ; ! k 2 )) and previous ( ^ ® i ¡ 1 ( ! k 1 ; ! k 2 )) estimates is smaller than a preselected threshold (e.g., ² = 10 ¡ 2 ): ¯ ¯ ¯ P K 1 ¡ 1 k 1 =0 P K 2 ¡ 1 k 2 =0 j ^ ® i ( ! k 1 ; ! k 2 ) j 2 ¡ P K 1 ¡ 1 k 1 =0 P K 2 ¡ 1 k 2 =0 j ^ ® i ¡ 1 ( ! k 1 ; ! k 2 ) j 2 ¯ ¯ ¯ P K 1 ¡ 1 k 1 =0 P K 2 ¡ 1 k 2 =0 j ^ ® i ¡ 1 ( ! k 1 ; ! k 2 ) j 2 · ²: (5.50) 5.5.1 Con v ergence Sp eed In our ¯rst example, w e study the con v ergence prop erties of the MAPES algorithms. W e use a 1-D example for simplicit y . (Note that a similar 1-D example

PAGE 79

67 w as considered in Chapter 3 but without the MAPES-CM algorithm whic h has b een in tro duced in this c hapter.) The true sp ectrum of the sim ulated signal is sho wn in Figure 3.1(a), where w e ha v e four sp ectral lines lo cated at f 1 = 0 : 05 Hz, f 2 = 0 : 065 Hz, f 3 = 0 : 26 Hz, and f 4 = 0 : 28 Hz with complex amplitudes ® 1 = ® 2 = ® 3 = 1 and ® 4 = 0 : 5. Besides these sp ectral lines, Figure 3.1(a) also sho ws a con tin uous sp ectral comp onen t cen tered at 0.18 Hz with a width of 0.015 Hz and a constan t mo dulus of 0.25. The data sequence has N 1 = 128 ( N 2 = 1) samples out of whic h 51 (40%) samples are missing; the lo cations of the missing samples are c hosen arbitrarily . The data are corrupted b y a zero-mean circularly symmetric complex white Gaussian noise with standard deviation 0.1. In Figure 5.1(a), the APES algorithm is applied to the complete data and the resulting sp ectrum is sho wn. The APES sp ectrum will b e used later as a reference for comparison purp oses. The WFFT sp ectrum for the incomplete data is sho wn in Figure 5.1(b). As exp ected, the WFFT sp ectrum has p o or resolution and high sidelob es, and it underestimates the true sp ectrum. Note that the WFFT sp ectrum will b e used as initial estimate for b oth the GAPES and MAPES algorithms in this example. Figure 5.1(c) sho ws the GAPES sp ectrum, whic h also underestimates the sin usoidal comp onen ts and giv es some artifacts, due to the p o or initial estimate of H ( ! 1 ; ! 2 ). The MAPES-CM sp ectrum is plotted in Figure 5.1(d). Figures. 5.1(e) and 5.1(f ) sho w the MAPES-EM1 and MAPES-EM2 sp ectra, resp ectiv ely . All the MAPES algorithms p erform quite w ell and their missing-data sp ectral estimates are similar to the high-resolution complete-data APES sp ectrum in Figure 5.1(a). The MAPES and GAPES sp ectral estimates at di®eren t iterations are plotted in Figures. 5.2(a) 5.2(d). Among the MAPES algorithms, MAPES-CM is the fastest to con v erge, after 3 iterations. MAPES-EM con v erges more slo wly , with MAPESEM1 con v erging after 11 and MAPES-EM2 after 9 iterations. Due to the relativ ely

PAGE 80

68 p o or initial ¯lter-bank H ( ! 1 ; ! 2 ) used in this example, the GAPES algorithm p erforms relativ ely p o orly and con v erges relativ ely slo wly , after 10 iterations. In the follo wing examples where the gapp ed-data initialization step in App endix 5.7 can b e applied, GAPES usually con v erges faster, within a few iterations. F or illustration purp oses, in Figure 5.2 w e only plot the ¯rst four iterations of eac h metho d and note that the con v ergence b eha vior of eac h algorithm do es not c hange signi¯can tly in the remaining iterations. 5.5.2 P erformance Study In this example, w e illustrate the p erformance of the MAPES algorithms for 2-D sp ectral estimation. W e consider a 16 £ 16 data matrix consisting of three 2-D sin usoids (signals 1, 2, and 3) at normalized frequencies (4/16, 5/16), (6/16, 5/16), and (10/16, 9/16) and with complex amplitudes equal to 1, 0.7, and 2, resp ectiv ely , em b edded in zero-mean circularly symmetric complex Gaussian white noise with standard deviation 0.1. All the samples in ro ws 4, 8, 11, 14, and in columns 3, 6, 7, 11, 12, 14 are missing, whic h amoun ts to o v er 50% of the total n um b er of data samples. The true amplitude sp ectrum is plotted in Figure 5.3(a) with the estimated amplitude v alues giv en next to eac h sin usoid. Eac h sp ectrum is obtained on a 64 £ 64 grid. The WFFT sp ectrum of the complete data is sho wn in Figure 5.3(b) along with the estimated magnitudes of the sin usoids at their true lo cations. The t w o closely spaced sin usoids are smeared together. Figure 5.3(c) sho ws the APES sp ectrum, also constructed from the complete data, whic h has w ell separated sp ectral p eaks and accurately estimated magnitudes. The data missing pattern is displa y ed in Figure 5.3(d). The WFFT sp ectrum for the missing data case is sho wn in Figure 5.3(e), whic h underestimates the sin usoids and con tains strong artifacts due to the zeros assumed for the missing samples. Figure 5.3(f ) sho ws the sp ectrum estimated b y GAPES, with an initial ¯lter of size 2 £ 2. In Figures. 5.3(g), 5.3(h), and 5.3(i), w e sho w the sp ectra estimated b y MAPES-EM1, MAPES-EM2, and MAPES-CM,

PAGE 81

69 resp ectiv ely . The GAPES algorithm estimates the signal magnitudes m uc h more accurately than the WFFT, but not as w ell as the MAPES algorithms. All MAPES algorithms p erform w ell giving accurate sp ectral estimates and clearly separated sp ectral p eaks. No w w e compare the computational complexities of these algorithms. The n um b ers of °oating p oin t op erations (Flops) needed for di®eren t algorithms are giv en in T able 5.1. The WFFT algorithm is the most e±cien t, and GAPES is more e±cien t than the MAPES algorithms. Among MAPES, MAPES-CM is 20 times faster than MAPES-EM2, whic h is ab out 17 times faster than MAPES-EM1. T able 5.1: Computational complexities of the WFFT, GAPES, MAPES-EM, and MAPES-CM sp ectral estimators. WFFT GAPES MAPES-EM1 MAPES-EM2 MAPES-CM Flops 3 £ 10 5 3 £ 10 9 1 £ 10 14 6 £ 10 12 3 £ 10 11 The results displa y ed so far w ere for one t ypical realization of the data. Using 100 Mon te Carlo sim ulations (b y v arying the realizations of the noise), w e obtain the ro ot mean-squared errors (RMSE's) of the magnitude and phase estimates of the three sin usoids at their true lo cations for WFFT, GAPES, and MAPES. The RMSE's of the magnitude and phase estimates are listed in T ables 5.2 and 5.3, resp ectiv ely . MAPES-EM1 giv es the b est accuracy follo w ed closely b y MAPES-EM2 and MAPESCM. GAPES is sligh tly less accurate, whereas the accuracy corresp onding to the WFFT is m uc h lo w er. 5.5.3 Syn thetic Ap erture Radar (SAR) Imaging Applications In the follo wing t w o examples, w e illustrate the applications of the missingdata algorithms to SAR imaging using incomplete phase-history data. 2-D high resolution phase-history data of a Slicy ob ject at 0 ± azim uth angle w ere generated b y using XP A TCH [2], a high frequency electromagnetic scattering

PAGE 82

70 T able 5.2: RMSE's of the magnitude estimates obtained via the WFFT, GAPES, MAPES-EM, and MAPES-CM sp ectral estimators. WFFT GAPES MAPES-EM1 MAPES-EM2 MAPES-CM Signal 1 0.577 0.021 0.010 0.013 0.014 Signal 2 0.372 0.053 0.009 0.014 0.014 Signal 3 1.156 0.011 0.008 0.011 0.015 T able 5.3: RMSE's of the phase (radian) estimates obtained via the WFFT, GAPES, MAPES-EM, and MAPES-CM sp ectral estimators. WFFT GAPES MAPES-EM1 MAPES-EM2 MAPES-CM Signal 1 0.148 0.010 0.009 0.009 0.009 Signal 2 0.145 0.013 0.012 0.015 0.016 Signal 3 0.057 0.006 0.004 0.004 0.006 prediction co de for complex 3-D ob jects. A photo of the Slicy ob ject tak en at 45 ± azim uth angle is sho wn in Figure 5.4(a). The original data matrix has a size of 288 £ 288 with a resolution of 0.043 meters in b oth range and cross-range. Figure 5.4(b) sho ws the mo dulus of the 2-D WFFT image of the original data. Here, w e consider only a 32 £ 32 cen ter blo c k of the phase-history data. The APES image of the complete 32 £ 32 data is sho wn in Figure 5.4(c). The data missing pattern is displa y ed in Figure 5.4(d), where the samples in ro ws 4, 14, 15, 22, 27, 28 and in columns 8, 9, 10, 20, 21, 22, 26, 27 are missing (p ossibly due to b oth angular div ersit y and strong in terferences), resulting in 40% missing data. Figure 5.4(e) sho ws the WFFT image, whic h has lo w resolution, high sidelob es, and smeared features. By using an initial ¯lter matrix of size 6 £ 6, the GAPES image is sho wn in Figure 5.4(f ), where strong artifacts around the dihedrals are readily observ ed. The image reconstructed via MAPES-CM is sho wn in Figure 5.4(g), whic h is quite similar to the complete-data image in Figure 5.4(c). Based on the data in terp olated via MAPES-CM, w e apply the recen tly dev elop ed rank-de¯cien t robust Cap on ¯lter-bank (R CF) sp ectral estimator

PAGE 83

71 [78 ] with a 20 £ 20 ¯lter-bank and a spherical steering v ector uncertain t y set with unit radius. The resulting image is sho wn in Figure 5.4(h), whic h exhibits no sidelob e problem and retains all imp ortan t features of Slicy . Compared with Figure 5.4(b), w e note that although the data size w as reduced from 288 £ 288 to 32 £ 32 and from the reduced data matrix 40% of the samples w ere omitted, w e can still obtain an image similar to that obtained b y the WFFT applied to the original high-resolution data. (Note that w e cannot get a similar high resolution image with all the w ell separated features and without sidelob e problem b y simply thresholding Figure 5.4(f ), 5.4(g), or ev en 5.4(c).) Next, w e consider a 32 £ 32 data matrix from the \Bac kho e Data Dome, V ersion 1.0.", whic h consists of sim ulated wideband (7-13 GHz), full p olarization, complex bac kscatter data from a bac kho e v ehicle in free space. The 3-D CAD mo del of the bac kho e v ehicle is sho wn in Figure 5.5(a). The bac kscattered data has b een generated o v er a full upp er 2 ¼ steradian viewing hemisphere, whic h is also illustrated in Figure 5.5(a). Figure 5.5(b) sho ws the 2-D in v erse FFT image (after the P olar-toCartesian in terp olation) formed at 0 ± elev ation angle with 6 GHz bandwidth and a 110 ± azim uth cut ranging from ¡ 10 ± to 100 ± . A t 0 ± elev ation, the reduced 32 £ 32 data matrix are collected from a 2 ± azim uth cut cen tered around 90 ± azim uth, co v ering a 0.3 GHz bandwidth cen tered around 10 GHz. Figure 5.5(c) sho ws the WFFT image of the complete data matrix with smeared features. Figure 5.5(d) sho ws the APES image of the complete data with a 16 £ 16 2-D ¯lter. Some smeared features in Figure 5.5(c) are clearly observ ed here, suc h as the one lo cated at ro w 26 and column 20. Figure 5.5(e) illustrates the data missing pattern, where the samples in ro w 5, 13, 14, 21, 22, 27 and in columns 8, 9, 10, 18, 19, 20, 26, 27 are missing. Figures 5.5(f ) and 5.5(g) sho w the WFFT and MAPES-CM images of the missing-data matrix. It can b e observ ed that despite the missing samples, MAPES-CM can still giv e a sp ectral estimate whic h has all the features sho wn in Figure 5.5(d).

PAGE 84

72 0 0.1 0.2 0.3 0.4 0.5 0 0.2 0.4 0.6 0.8 1 1.2 Frequency (Hz)Modulus of Complex Amplitude 0 0.1 0.2 0.3 0.4 0.5 0 0.2 0.4 0.6 0.8 1 1.2 Frequency (Hz)Modulus of Complex Amplitude (a) (b) 0 0.1 0.2 0.3 0.4 0.5 0 0.2 0.4 0.6 0.8 1 1.2 Frequency (Hz)Modulus of Complex Amplitude 0 0.1 0.2 0.3 0.4 0.5 0 0.2 0.4 0.6 0.8 1 1.2 Frequency (Hz)Modulus of Complex Amplitude (c) (d) 0 0.1 0.2 0.3 0.4 0.5 0 0.2 0.4 0.6 0.8 1 1.2 Frequency (Hz)Modulus of Complex Amplitude 0 0.1 0.2 0.3 0.4 0.5 0 0.2 0.4 0.6 0.8 1 1.2 Frequency (Hz)Modulus of Complex Amplitude (e) (f ) Figure 5.1: Mo dulus of the missing-data sp ectral estimates ( N 1 = 128, N 2 = 1, 51 (40%) missing samples): (a) Complete-data APES, (b) WFFT, (c) GAPES with M 1 = 64 and M 2 = 1, (d) MAPES-CM with M 1 = 64 and M 2 = 1, (e) MAPES-EM1 with M 1 = 64 and M 2 = 1, and (f ) MAPES-EM2 with M 1 = 64 and M 2 = 1.

PAGE 85

73 0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 Modulus i=0 0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 Modulus i=1 0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 Modulus i=2 0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 Modulus i=3 0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 Modulus i=4 Frequency (Hz) 0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 Modulus i=0 0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 Modulus i=1 0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 Modulus i=2 0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 Modulus i=3 0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 Modulus i=4 Frequency (Hz) (a) (b) 0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 Modulus i=0 0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 Modulus i=1 0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 Modulus i=2 0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 Modulus i=3 0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 Modulus i=4 Frequency (Hz) 0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 Modulus i=0 0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 Modulus i=1 0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 Modulus i=2 0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 Modulus i=3 0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 Modulus i=4 Frequency (Hz) (c) (d) Figure 5.2: Mo dulus of the missing-data sp ectral estimates obtained at di®eren t iterations ( N 1 = 128, N 2 = 1, 51 (40%) missing samples): (a) MAPES-CM, (b) MAPES-EM1, (c) MAPES-EM2 , and (d) GAPES.

PAGE 86

74 1.0000 0.7000 2.0000 0 0.2 0.4 0.6 0.8 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.99808 0.68569 2.0014 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 1.005 0.70055 2.0017 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 (a) (b) (c) 2 4 6 8 10 12 14 16 2 4 6 8 10 12 14 16 0.42636 0.32908 0.84218 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 0.9768 0.64683 1.9866 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 (d) (e) (f ) 1.0069 0.70073 1.9916 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 1.0136 0.69643 1.9936 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 0.99644 0.70751 1.9853 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 (g) (h) (i) Figure 5.3: Mo dulus of the 2-D sp ectra: (a) T rue sp ectrum, (b) 2-D complete-data WFFT, (c) 2-D complete-data APES with M 1 = M 2 = 8, (d) 2-D data missing pattern, the blac k strip es indicate missing samples, (e) 2-D WFFT, (f ) 2-D GAPES with M 1 = M 2 = 8, (g) 2-D MAPES-EM1 with M 1 = M 2 = 8, (h) 2-D MAPES-EM2 with M 1 = M 2 = 8, and (i) 2-D MAPES-CM with M 1 = M 2 = 8.

PAGE 87

75 (a) (b) 5 10 15 20 25 30 5 10 15 20 25 30 (c) (d) (e) (f ) (g) Figure 5.4: Mo dulus of the SAR images of the Slicy ob ject obtained from a 32 £ 32 data matrix with missing samples: (a) Photograph of the ob ject (tak en at 45 ± azim uth angle), (b) 2-D WFFT for a 288 £ 288 (not 32 £ 32) data matrix, (c) 2-D completedata APES with M 1 = M 2 = 16, (d) 2-D data missing pattern, the blac k strip es indicate missing samples, (e) 2-D WFFT, (f ) 2-D GAPES with M 1 = M 2 = 16, and (g) 2-D MAPES-CM with M 1 = M 2 = 16.

PAGE 88

76 (a) 100 200 300 400 500 600 700 800 900 1000 100 200 300 400 500 600 700 800 900 1000 10 20 30 40 50 60 10 20 30 40 50 60 10 20 30 40 50 60 10 20 30 40 50 60 (b) (c) (d) 5 10 15 20 25 30 5 10 15 20 25 30 10 20 30 40 50 60 10 20 30 40 50 60 10 20 30 40 50 60 10 20 30 40 50 60 (e) (f ) (g) Figure 5.5: Mo dulus of the SAR images of the Bac kho e data obtained from a 32 £ 32 data matrix with missing samples: (a) 3-D CAD mo del of the target and the 3-D K-space data dome, (b) mo dulus of the 2-D FFT image formed at 0 ± elev ation with 6 GHz bandwidth and a 110 ± azim uth cut (not 32 £ 32), (c) 2-D complete-data WFFT, (d) 2-D complete-data APES with a 2-D ¯lter of size 16 £ 16, (e) 2-D data missing pattern, the blac k strip es indicate missing samples, (f ) 2-D WFFT, and (g) 2-D MAPES-CM with a 2-D ¯lter of size 16 £ 16.

PAGE 89

77 5.6 Conclusions W e ha v e considered 2-D nonparametric complex sp ectral estimation (with its 1-D coun terpart as a sp ecial case) for data matrices with missing samples o ccurring in arbitrary patterns. The previously prop osed MAPES-EM algorithms ha v e b een extended to the 2-D case. W e ha v e also dev elop ed another missing-data algorithm, referred to as MAPES-CM, b y maximizing an ML ¯tting criterion iterativ ely via a cyclic maximization (CM) algorithm. W e ha v e compared MAPES-EM with MAPESCM and ha v e sho wn that MAPES-CM allo ws signi¯can t computational sa vings compared with MAPES-EM, whic h is esp ecially desirable for 2-D applications, and also that MAPES-CM con v erges faster than MAPES-EM. Numerical examples ha v e b een pro vided to demonstrate the p erformance of the MAPES algorithms. 5.7 App endix: 2-D GAPES In this app endix, w e sho w the extension of the GAPES algorithm to t w odimensional (2-D) data matrices. Let G b e the set of sample indices ( n 1 ; n 2 ) for whic h the data samples are a v ailable, and U b e the set of sample indices ( n 1 ; n 2 ) for whic h the data samples are missing. The set of a v ailable samples f y n 1 ;n 2 : ( n 1 ; n 2 ) 2 G g is denoted b y the g £ 1 v ector ° , whereas the set of missing samples f y n 1 ;n 2 : ( n 1 ; n 2 ) 2 U g is denoted b y the ( N 1 N 2 ¡ g ) £ 1 v ector ¹ . The problem of in terest is to estimate ® ( ! 1 ; ! 2 ) giv en ° . Assume w e consider a K 1 £ K 2 -p oin t DFT grid: ( ! k 1 ; ! k 2 ) = (2 ¼ k 1 =K 1 ; 2 ¼ k 2 = K 2 ) for k 1 = 0 ; :::; K 1 ¡ 1 and k 2 = 0 ; :::; K 2 ¡ 1 (with K 1 > N 1 and K 2 > N 2 ). The 2-D GAPES algorithm tries to solv e the follo wing minimization problem: min ¹ ; f ® ( ! k 1 ;! k 2 ) ; H ( ! k 1 ;! k 2 ) g P K 1 ¡ 1 k 1 =0 P K 2 ¡ 1 k 2 =0 k v ec( H ( ! k 1 ; ! k 2 ) ? Y ) ¡ ® ( ! k 1 ; ! k 2 ) a L 1 ;L 2 ( ! k 1 ; ! k 2 ) k 2 sub ject to v ec H ( H ( ! k 1 ; ! k 2 )) a M 1 ;M 2 ( ! k 1 ; ! k 2 ) = 1 ; (5.51) via cyclic optimization.

PAGE 90

78 F or the initialization step, w e obtain the initial APES estimates of H ( ! 1 ; ! 2 ) and ® ( ! 1 ; ! 2 ) from the a v ailable data ° in the follo wing w a y . Let S b e the set of snapshot indices ( l 1 ; l 2 ) suc h that the elemen ts of the corresp onding initial data snapshot indices f ( l 1 ; l 2 ) ; :::; ( l 1 ; l 2 + M 0 2 ¡ 1) ; :::; ( l + M 0 1 ¡ 1 ; l 2 ) ; :::; ( l 1 + M 0 1 ¡ 1 ; l 2 + M 0 2 ¡ 1) g 2 G . De¯ne the set of M 0 1 M 0 2 £ 1 v ectors f ¹ y l 1 ;l 2 : ( l 1 ; l 2 ) 2 S g whic h con tain a v ailable data samples only and let jS j b e the n um b er of v ectors in S . F urthermore, de¯ne the initial sample co v ariance matrix ^ R = 1 jS j X ( l 1 ;l 2 ) 2S ¹ y l 1 ;l 2 ¹ y H l 1 ;l 2 : (5.52) The size of the initial ¯lter matrix M 0 1 £ M 0 2 m ust b e c hosen suc h that the ^ R calculated in (5.52) has full rank. Similarly , the initial F ourier transform of the data snapshots is giv en b y ¹ g ( ! 1 ; ! 2 ) = 1 jS j X ( l 1 ;l 2 ) 2S ¹ y l 1 ;l 2 e ¡ j ( ! 1 l 1 + ! 2 l 2 ) : (5.53) So the initial estimates of H ( ! 1 ; ! 2 ) and ® ( ! 1 ; ! 2 ) can b e calculated b y (4.13) (4.15) but b y using the ^ R and ¹ g ( ! 1 ; ! 2 ) giv en ab o v e. Next, w e in tro duce some additional notation that will b e used later for the step of in terp olating the missing samples. Let the L 1 L 2 £ ( L 2 N 1 ¡ M 1 + 1) matrix T b e de¯ned b y T = 2 6 6 6 4 I L 1 0 L 1 ;M 1 ¡ 1 I L 1 0 L 1 ;M 1 ¡ 1 . . . I L 1 3 7 7 7 5 : (5.54) Hereafter, 0 K 1 ;K 2 denotes a K 1 £ K 2 matrix of zeros only and I K stands for the K £ K iden tit y matrix. F urthermore, let G b e the follo wing ( L 2 N 1 ¡ M 1 + 1) £ N 1 N 2 T o eplitz matrix G ( ! 1 ; ! 2 ) = 2 6 6 6 4 h H 1 0 1 ;L 1 ¡ 1 h H 2 0 1 ;L 1 ¡ 1 ¢ ¢ ¢ h H M 2 0 ¢ ¢ ¢ 0 0 h H 1 0 1 ;L 1 ¡ 1 h H 2 0 1 ;L 1 ¡ 1 ¢ ¢ ¢ h H M 2 ¢ ¢ ¢ 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 ¢ ¢ ¢ 0 h H 1 0 1 ;L 1 ¡ 1 h H 2 0 1 ;L 1 ¡ 1 ¢ ¢ ¢ h H M 2 3 7 7 7 5 ; (5.55)

PAGE 91

79 where f h m 2 g M 2 m 2 =1 are the corresp onding columns of H ( ! 1 ; ! 2 ). With these de¯nitions, w e ha v e v ec ( ¹ X ) = v ec( H ( ! 1 ; ! 2 ) ? Y ) = TG v ec ( Y ) : (5.56) By making use of (5.56), the estimate of ¹ based on the initial estimates ^ ® ( ! 1 ; ! 2 ) and ^ H ( ! 1 ; ! 2 ) is giv en b y the solution to the follo wing problem min ¹ ° ° ° ° ° ° ° 2 6 4 TG ( ! 0 ; ! 0 ) . . . TG ( ! K 1 ¡ 1 ; ! K 2 ¡ 1 ) 3 7 5 v ec ( Y ) ¡ 2 6 4 ^ ® ( ! 0 ; ! 0 ) a L 1 ;L 2 ( ! 0 ; ! 0 ) . . . ^ ® ( ! K 1 ¡ 1 ; ! K 2 ¡ 1 ) a L 1 ;L 2 ( ! K 1 ¡ 1 ; ! K 2 ¡ 1 ) 3 7 5 ° ° ° ° ° ° ° 2 : (5.57) T o solv e (5.57), let the matrices G ° ( ! k 1 ; ! k 2 ) and G ¹ ( ! k 1 ; ! k 2 ) b e de¯ned implicitly b y the follo wing equalit y: G ( ! k 1 ; ! k 2 )v ec( Y ) = G ° ( ! k 1 ; ! k 2 ) ° + G ¹ ( ! k 1 ; ! k 2 ) ¹ ; (5.58) where ° and ¹ are the v ectors con taining the a v ailable samples and missing samples, resp ectiv ely . In other w ords, G ° ( ! k 1 ; ! k 2 ) and G ¹ ( ! k 1 ; ! k 2 ) con tain the columns of G ( ! k 1 ; ! k 2 ) that corresp ond to the indices in G and U , resp ectiv ely . By in tro ducing the follo wing matrices: ~ G ° 4 = 2 6 4 TG ° ( ! 0 ; ! 0 ) . . . TG ° ( ! K 1 ¡ 1 ; ! K 2 ¡ 1 ) 3 7 5 (5.59) and ~ G ¹ 4 = 2 6 4 TG ¹ ( ! 0 ; ! 0 ) . . . TG ¹ ( ! K 1 ¡ 1 ; ! K 2 ¡ 1 ) 3 7 5 ; (5.60) the criterion (5.57) can then b e written as min ¹ ° ° ° ~ G ¹ ¹ + ~ G ° ° ¡ ~ ® ° ° ° 2 ; (5.61) where ~ ® 4 = 2 6 4 ^ ® ( ! 0 ; ! 0 ) a L 1 ;L 2 ( ! 0 ; ! 0 ) . . . ^ ® ( ! K 1 ¡ 1 ; ! K 2 ¡ 1 ) a L 1 ;L 2 ( ! K 1 ¡ 1 ; ! K 2 ¡ 1 ) 3 7 5 (5.62)

PAGE 92

80 The closed-form solution of the quadratic problem (5.61) is easily obtained as ^ ¹ = ³ ~ G H ¹ ~ G ¹ ´ ¡ 1 ~ G H ¹ ³ ~ ® ¡ ~ G ° ° ´ : (5.63) A step-b y-step summary of 2-D GAPES follo ws. Step 0: Obtain an initial estimate of f ® ( ! 1 ; ! 2 ) ; h ( ! 1 ; ! 2 ) g . Step 1: Use the most recen t estimate of f ® ( ! 1 ; ! 2 ) ; h ( ! 1 ; ! 2 ) g in (5.51) to estimate ¹ b y minimizing the so-obtained cost function, whose solution is giv en b y (5.63). Step 2: Use the latest estimate of ¹ to ¯ll in the missing data samples, and estimate f ® ( ! 1 ; ! 2 ) ; H ( ! 1 ; ! 2 ) g K 1 ¡ 1 ;K 2 ¡ 1 k 1 =0 ;k 2 =0 b y minimizing the cost function in (5.51) based on the in terp olated data. (This step is equiv alen t to applying 2-D APES to the complete data.) Step 3: Rep eat Steps 1 to 2 un til practical con v ergence.

PAGE 93

CHAPTER 6 RANK-DEFICIENT R CF SPECTRAL ESTIMA TOR 6.1 In tro duction All the missing-data algorithms dev elop ed in the previous c hapters estimate the missing data samples. This indicates that b etter sp ectral estimation p erformance, e.g. higher resolution than APES, could b e ac hiev ed based on the in terp olated (complete) data. It is w ell kno wn that Cap on sp ectral estimator has b etter resolution (narro w er mainlob e) compared with APES. In this c hapter, w e dev elop a Cap on-based sp ectral estimator whic h can ac hiev e higher resolution than the MAPES algorithms. The adaptiv e FIR ¯lter-bank used in the Cap on sp ectral estimator is obtained via the Cap on b eamformer [71 , 77], whic h w as originally prop osed for adaptiv e b eamforming [10 , 28 ]. (W e \b orro w" the terminology \b eamforming" for sp ectral estimation in order to quic kly imp ort results from arra y pro cessing.) The data-dep enden t Cap on b eamformer adaptiv ely selects the w eigh t v ector to minimize the arra y output p o w er sub ject to the linear constrain t that the signal-of-in terest (SOI) do es not su®er from an y distortion [10 , 28 ]. The Cap on b eamformer has b etter resolution and m uc h b etter in terference rejection capabilit y than the data-indep enden t b eamformer, pro vided that the arra y steering v ector corresp onding to the SOI is accurately kno wn and that the n um b er of snapshots is relativ ely large. Ho w ev er, for complex sp ectral estimation, the ¯lter length is often c hosen to b e quite large in order to ac hiev e high resolution; hence the n um b er of \snapshots" is usually small. Whenev er this happ ens, the Cap on b eamformer ma y suppress the SOI as if it w ere an in terference, whic h results in a signi¯can tly underestimated SOI p o w er. This is, in fact, the reason wh y the Cap on sp ectral estimates are generally biased do wn w ard. 81

PAGE 94

82 Man y approac hes ha v e b een prop osed during the past three decades to impro v e the robustness of the Cap on b eamformer. See, for example, [77 , 38 ] and the references therein for the extensiv e literature on robust adaptiv e b eamforming. In [38 , 75 ], w e obtained a Robust Cap on Beamformer (R CB) via a co v ariance matrix ¯tting approac h [43 ]. The computational complexit y of R CB is comparable to that of the standard Cap on b eamformer. In this c hapter, w e will sho w that R CB allo ws the sample co v ariance matrix to b e rank-de¯cien t and hence it can b e used to design an FIR ¯lter-bank to ac hiev e ev en higher resolution than the standard Cap on approac h for complex sp ectral estimation. Using rank-de¯cien t sample co v ariance matrices for sp ectral estimation w as ¯rst considered b y Benitz (see, e.g., [3 ] and the references therein). In particular, he referred to using suc h sp ectral estimation metho ds for SAR image formation as high-de¯nition imaging (HDI). It has b een sho wn in [47 ] that HDI can b e used to signi¯can tly impro v e the automatic target recognition (A TR) p erformance of a mo dern SAR, whic h demonstrates the imp ortance of sp ectral estimation based on rank-de¯cien t sample co v ariance matrices. In this c hapter, w e consider nonparametric complex sp ectral estimation using an adaptiv e ¯ltering based approac h where the FIR ¯lter-bank is obtained via a rankde¯cien t R CB. W e deriv e the rank-de¯cien t robust Cap on ¯lter-bank (R CF) sp ectral estimator in detail. W e sho w that b y allo wing the sample co v ariance matrix to b e rank-de¯cien t, w e can ac hiev e m uc h higher resolution than existing approac hes, whic h is useful in man y applications including radar target detection and feature extraction. Numerical examples are pro vided to demonstrate the p erformance of the new approac h as compared to data-adaptiv e and data-indep enden t FIR ¯ltering based sp ectral estimation metho ds. The remainder of this c hapter is organized as follo ws. The problem of in terest is form ulated in Section 6.2. The deriv ation of the rank-de¯cien t R CF complex

PAGE 95

83 sp ectral estimator is giv en in Section 6.3. Both one-dimensional (1-D) and t w odimensional (2-D) n umerical examples are pro vided in Section 6.4 to illustrate the p erformance of the rank-de¯cien t R CF for sp ectral estimation and SAR imaging as compared with other approac hes. Finally , Section 6.5 concludes this c hapter. 6.2 Problem F orm ulation and Preliminaries Consider the problem of estimating the amplitude sp ectrum of a complexv alued discrete-time 1-D signal f y n g N ¡ 1 n =0 . (Extension to 2-D data, whic h will b e used in one of the n umerical examples in Section 6.4, is giv en in App endix 6.6.1.) F or a frequency ! of in terest, the signal y n is mo deled as y n = ® ( ! ) e j ! n + e n ( ! ) ; n = 0 ; :::; N ¡ 1 ; ! 2 [0 ; 2 ¼ ) ; (6.1) where ® ( ! ) denotes the complex amplitude of the sin usoidal comp onen t with frequency ! , and e n ( ! ) denotes the residual term (assumed zero-mean) at frequency ! , whic h includes the unmo deled noise and in terference from frequencies other than ! . The problem of in terest is to estimate ® ( ! ) from f y n g N ¡ 1 n =0 for an y giv en frequency ! . The ¯lter-bank approac hes address the aforemen tioned sp ectral estimation problem b y passing the data f y n g through a bank of FIR bandpass ¯lters with v arying cen ter frequency ! , and then obtaining the amplitude sp ectrum estimate ^ ® ( ! ) for ! 2 [0 ; 2 ¼ ) from the ¯ltered data. W e denote an M -tap FIR bandpass ¯lter b y h ( ! ) = [ h 0 h 1 ¢ ¢ ¢ h M ¡ 1 ] T ; (6.2) where ( ¢ ) T denotes the transp ose. Let the forw ard data v ectors ¹ y l = [ y l y l +1 ¢ ¢ ¢ y l + M ¡ 1 ] T ; l = 0 ; :::; L ¡ 1 (6.3) b e the o v erlapping M £ 1 sub v ectors constructed from the data v ector ¹ y = [ y 0 y 1 ¢ ¢ ¢ y N ¡ 1 ] T ; (6.4)

PAGE 96

84 where L = N ¡ M + 1. Then, according to the data mo del in (6.1), the forw ard data v ectors can b e written as ¹ y l = ® ( ! ) a ( ! ) ¢ e j ! l + ¹ e l ( ! ) ; (6.5) where a ( ! ) is an M £ 1 v ector giv en b y a ( ! ) = [1 e j ! ¢ ¢ ¢ e j ! ( M ¡ 1) ] T (6.6) and ¹ e l ( ! ) = [ e l ( ! ) e l +1 ( ! ) ¢ ¢ ¢ e l + M ¡ 1 ( ! )] T . Hence the output samples obtained b y passing ¹ y l through the FIR ¯lter h ( ! ) can b e written as h H ( ! ) ¹ y l = ® ( ! )[ h H ( ! ) a ( ! )] e j ! l + ¹ w l ( ! ) ; (6.7) where ( ¢ ) H denotes the conjugate transp ose and ¹ w l ( ! ) = h H ( ! ) ¹ e l ( ! ) denotes the residue term at the ¯lter output. F or an undistorted sp ectral estimate, w e require that h H ( ! ) a ( ! ) = 1 : (6.8) F rom the output of the FIR ¯lter h H ( ! ) ¹ y l = ® ( ! ) e j ! l + ¹ w l ( ! ) ; (6.9) w e can obtain the least-squares estimate of ® ( ! ) as ^ ¹ ® ( ! ) = h H ( ! ) ¹ g ( ! ) ; (6.10) where ¹ g ( ! ) is the normalized F ourier transform of the forw ard data v ectors ¹ g ( ! ) = 1 L L ¡ 1 X l =0 ¹ y l e ¡ j ! l : (6.11) Since a com bined forw ard-bac kw ard approac h usually yields b etter sp ectral estimates than the forw ard-only approac h, w e also consider the bac kw ard data v ectors ~ y l = [ y ¤ N ¡ l ¡ 1 y ¤ N ¡ l ¡ 2 ¢ ¢ ¢ y ¤ N ¡ l ¡ M ] T ; l = 0 ; :::; L ¡ 1 ; (6.12)

PAGE 97

85 where ( ¢ ) ¤ denotes the complex conjugate. Note that f ~ y l g L ¡ 1 l =0 are the o v erlapping M £ 1 sub v ectors constructed from the data v ector ~ y = [ y ¤ N ¡ 1 y ¤ N ¡ 2 ¢ ¢ ¢ y ¤ 0 ] T : (6.13) Similarly to ¹ y l , ~ y l can b e written as ~ y l = ® ¤ ( ! ) e ¡ j ( N ¡ 1) ! a ( ! ) ¢ e j ! l + ~ e l ( ! ) ; (6.14) where ~ e l ( ! ) = [ e ¤ N ¡ l ¡ 1 ( ! ) e ¤ N ¡ l ¡ 2 ( ! ) ¢ ¢ ¢ e ¤ N ¡ l ¡ M ( ! )] T . P assing ~ y l through the FIR ¯lter h ( ! ) yields the follo wing output h H ( ! ) ~ y l = e ¡ j ( N ¡ 1) ! ¢ ® ¤ ( ! ) e j ! l + ~ w l ( ! ) ; (6.15) where ~ w l ( ! ) = h H ( ! ) ~ e l ( ! ) denotes the residue term at the ¯lter output. F rom the ab o v e FIR ¯lter output, w e can obtain another least-squares estimate of ® ( ! ): ^ ~ ® = e ¡ j ( N ¡ 1) ! ~ g H ( ! ) h ( ! ) ; (6.16) where ~ g ( ! ) is the normalized F ourier transform of the bac kw ard data v ectors: ~ g ( ! ) = 1 L L ¡ 1 X l =0 ~ y l e ¡ j ! l : (6.17) Av eraging the t w o least-squares estimates, ^ ¹ ® ( ! ) and ^ ~ ® ( ! ), giv es the forw ardbac kw ard estimate of ® ( ! ): ^ ® ( ! ) = 1 2 [ h H ( ! ) ¹ g ( ! ) + e ¡ j ( N ¡ 1) ! ~ g H ( ! ) h ( ! )] : (6.18) The forw ard-bac kw ard approac h is used in all of the adaptiv e ¯ltering based sp ectral estimators in the sections to follo w, and also in the determination of h ( ! ) whic h is discussed next. 6.3 Rank-De¯cien t Robust Cap on Filter-Bank Sp ectral Estimator W e deriv e the rank-de¯cien t R CF sp ectral estimator in a co v ariance matrix ¯tting framew ork b y assuming that the sample co v ariance matrix is singular, whic h

PAGE 98

86 happ ens, for example, when the FIR ¯lter length M is large. In particular, w e use the rank-de¯cien t R CB approac h to determine the data-dep enden t FIR ¯lter h ( ! ) from the sample co v ariance matrix. Hence, b esides in tro ducing a high-resolution sp ectral estimator, our deriv ations will also shed more ligh t on the prop erties of the R CB algorithm when the sample co v ariance matrix is singular. Note that, in principle, the same rank-de¯cien t approac h can b e used in arra y signal pro cessing whenev er the n um b er of a v ailable snapshots is less than the n um b er of sensors. 6.3.1 Rank-De¯cien t Sample Co v ariance Matrix It can b e observ ed from Section 6.2 that the forw ard and bac kw ard data v ectors are related b y ~ y l = J ¹ y ¤ L ¡ l ¡ 1 ; (6.19) where J denotes the exc hange matrix whose an tidiagonal elemen ts are ones and all the others are zeros. Similarly , for eac h frequency ! of in terest, w e ha v e ~ e l ( ! ) = J ¹ e ¤ L ¡ l ¡ 1 ( ! ) : (6.20) Let the T o eplitz noise co v ariance matrix Q ( ! ) b e de¯ned b y Q ( ! ) 4 = E £ ¹ e l ( ! ) ¹ e H l ( ! ) ¤ = E £ ~ e l ( ! ) ~ e H l ( ! ) ¤ ; (6.21) where the second equalit y follo ws from (6.20) and the T o eplitz structure of Q ( ! ). The co v ariance matrix of ¹ y l or, equiv alen tly , of ~ y l is giv en b y , R = j ® ( ! ) j 2 a ( ! ) a H ( ! ) + Q ( ! ) : (6.22) Let ^ R f and ^ R b denote the sample co v ariance matrices estimated from f ¹ y l g and f ~ y l g , resp ectiv ely , as follo ws: ^ R f = 1 L L ¡ 1 X l =0 ¹ y l ¹ y H l (6.23) ^ R b = 1 L L ¡ 1 X l =0 ~ y l ~ y H l : (6.24)

PAGE 99

87 Then the forw ard-bac kw ard estimate of the co v ariance matrix R is giv en b y ^ R = 1 2 ( ^ R f + ^ R b ) : (6.25) F rom (6.19), it is straigh tforw ard to sho w that ^ R b = J ^ R T f J ; (6.26) and hence, the ^ R in (6.25) is p ersymmetric. Compared with the non-p ersymmetric sample co v ariance matrix ^ R f estimated only from the forw ard data v ectors f ¹ y l g , the forw ard-bac kw ard ^ R is generally a b etter estimate of the true R . The dataadaptiv e FIR ¯ltering based sp ectral estimation metho ds w e consider herein obtain the data-dep enden t FIR ¯lter h ( ! ) from the ab o v e ^ R . Let ^ R b e the M £ M p ositiv e semi-de¯nite sample co v ariance matrix de¯ned in (6.25). Let K denote the rank of ^ R . With probabilit y one, K = 2 L , assuming that M > 2 L or, equiv alen tly , M > 2 3 ( N + 1). By c ho osing suc h a large M , w e hop e to ac hiev e high resolution for sp ectral estimation. Let ^ R = ^ S ^ ª ^ S H (6.27) denote the eigen-decomp osition of ^ R , where ^ S is an M £ K semi-unitary matrix with full column rank ( K < M ) and ^ ª is a K £ K p ositiv e de¯nite diagonal matrix whose diagonal elemen ts are the eigen v alues of ^ R . Next, w e deriv e the rank-de¯cien t R CF sp ectral estimator based on this singular ^ R . 6.3.2 Robust Cap on Filter-Bank (R CF) Approac h Owing to the small snapshot n um b er problem, the signal term in ^ R is not w ell describ ed b y j ® ( ! ) j 2 a ( ! ) a H ( ! ), but b y j ® ( ! ) j 2 ^ a ( ! ) ^ a H ( ! ) with ^ a ( ! ) b eing some v ector in the vicinit y of a ( ! ) and ^ a ( ! ) 6 = a ( ! ) [15 ]. (More generally , the smallsample errors of the co v ariance matrix can b e view ed as steering v ector errors in a corresp onding theoretical co v ariance matrix, and vice v ersa, see App endix 6.6.4;

PAGE 100

88 also see [72 ] for more details.) Consequen tly , if w e designed h ( ! ) b y means of the standard Cap on b eamformer: min h ( ! ) h H ( ! ) ^ R h ( ! ) sub ject to h H ( ! ) a ( ! ) = 1 ; (6.28) the Euclidean norm of h ( ! ), denoted k h ( ! ) k , w ould result rather large since ^ a ( ! ) is close to a ( ! ) and h ( ! ) passes a ( ! ) but attempts to suppress ^ a ( ! ). A large k h ( ! ) k indicates a large noise gain, whic h ma y sev erely degrade the estimation accuracy of ® ( ! ). It follo ws that w e should design h ( ! ) b y min h ( ! ) h H ( ! ) ^ R h ( ! ) sub ject to h H ( ! ) ^ a ( ! ) = ½; (6.29) where ½ is determined b y the constrain t h H ( ! ) a ( ! ) = 1; see (6.8) (note that ½ will b e close to 1 since ^ a ( ! ) is close to a ( ! )). In summary , w e design h ( ! ) based on ^ a ( ! ), instead of a ( ! ), to a v oid a large noise gain; furthermore, w e c ho ose ½ based on the constrain t h H ( ! ) a ( ! ) = 1 to get an un biased estimate of ® ( ! ) when w e use h ( ! ) in (6.18). The solution to (6.29) is deriv ed as follo ws. 6.3.2.1 Assuming that ^ a ( ! ) is giv en With · h ( ! ) = h ( ! ) =½ , (6.29) can b e rewritten as min · h ( ! ) · h H ( ! ) · R · h ( ! ) sub ject to · h H ( ! ) ^ a ( ! ) = 1 ; (6.30) where · R 4 = j ½ j 2 ^ R = ^ S · ª ^ S H and · ª 4 = j ½ j 2 ^ ª . Let ^ G denote the M £ ( M ¡ K ) matrix whose columns are the eigen v ectors of ^ R (or · R ) corresp onding to the zero eigen v alues. Hence ^ G spans the orthogonal complemen t of the subspace spanned b y ^ S . Let · h ( ! ) b e written as · h ( ! ) = ^ S · h 1 ( ! ) + ^ G · h 2 ( ! ) ; (6.31) where · h 1 ( ! ) = ^ S H · h ( ! ) and · h 2 ( ! ) = ^ G H · h ( ! ). ² Case 1: ^ a ( ! ) b elongs to the range space of · R .

PAGE 101

89 Let ^ a ( ! ) b e written as ^ a ( ! ) 4 = ^ S ³ for some non-zero v ector ³ . Using (6.31) in (6.30) for this case yields min · h 1 ( ! ) · h H 1 ( ! ) · ª · h 1 ( ! ) sub ject to · h H 1 ( ! ) ³ = 1 ; (6.32) whic h can b e readily solv ed as · h 1 ( ! ) = · ª ¡ 1 ³ ³ H · ª ¡ 1 ³ : (6.33) Since · h 2 ( ! ) is irrelev an t in this case, w e let · h 2 ( ! ) = 0 (6.34) to reduce the noise gain of · h ( ! ) ( k · h ( ! ) k 2 = k · h 1 ( ! ) k 2 + k · h 2 ( ! ) k 2 ). Then the FIR ¯lter has the form · h ( ! ) = ^ S · ª ¡ 1 ³ ³ H · ª ¡ 1 ³ = · R y ^ a ( ! ) ^ a H ( ! ) · R y ^ a ( ! ) ; (6.35) where · R y = ^ S · ª ¡ 1 ^ S H is the Mo ore-P enrose pseudo-in v erse of · R . Consequen tly , ^ h ( ! ) = ½ ¢ ^ R y ^ a ( ! ) ^ a H ( ! ) ^ R y ^ a ( ! ) : (6.36) Substituting (6.36) in to h H ( ! ) a ( ! ) = 1, w e ha v e ½ = ^ a H ( ! ) ^ R y ^ a ( ! ) a H ( ! ) ^ R y ^ a ( ! ) : (6.37) Hence ^ h ( ! ) is giv en b y ^ h ( ! ) = ^ R y ^ a ( ! ) a H ( ! ) ^ R y ^ a ( ! ) (6.38) and w e obtain the complex sp ectral estimator b y substituting (6.38) in to (6.18): ^ ® ( ! ) = 1 2 [ ^ h H ( ! ) ¹ g ( ! ) + e ¡ j ( N ¡ 1) ! ~ g H ( ! ) ^ h ( ! )] : (6.39) ² Case 2: ^ a ( ! ) do es not b elong to the range space of · R .

PAGE 102

90 Let ^ a ( ! ) b e written as ^ a ( ! ) = ^ S ³ + ^ G ¯ for some non-zero v ectors ³ and ¯ . No w (6.30) b ecomes min · h 1 ( ! ) ; · h 2 ( ! ) · h H 1 ( ! ) · ª · h 1 ( ! ) sub ject to · h H 1 ( ! ) ³ + · h H 2 ( ! ) ¯ = 1 ; (6.40) whic h admits a trivial solution: · h 1 ( ! ) = 0 (6.41) and (for example) · h 2 ( ! ) = ¯ k ¯ k 2 : (6.42) Consequen tly , since ¹ g ( ! ) and ~ g ( ! ) are linear transformations of f ¹ y l g and f ~ y l g , resp ectiv ely , where f ¹ y l g and f ~ y l g are in the range space of · R , w e ha v e · h H ( ! ) ¹ g ( ! ) = · h H ( ! ) ~ g ( ! ) = 0 ; (6.43) whic h giv es ^ ® ( ! ) = 0 : (6.44) Com bining the aforemen tioned t w o cases, the complex sp ectral estimate can b e written as ^ ® ( ! ) = ½ 1 2 [ ^ h H ( ! ) ¹ g ( ! ) + e ¡ j ( N ¡ 1) ! ~ g H ( ! ) ^ h ( ! )] ; ^ a ( ! ) 2 R ( ^ R ) 0 ; ^ a ( ! ) = 2 R ( ^ R ) ; (6.45) where ^ h ( ! ) is giv en in (6.38) and R ( ^ R ) denotes the range space of ^ R . W e remark that w e could ha v e obtained a signal p o w er estimate from (6.30) as follo ws: ^ · ¾ 2 ( ! ) 4 = j ^ ® ( ! ) j 2 = ^ h H ( ! ) ^ R ^ h ( ! ) = · h H ( ! ) · R · h ( ! ) = ( j ½ j 2 ^ a H ( ! ) ^ R y ^ a ( ! ) ; ^ a ( ! ) 2 R ( ^ R ) 0 ; ^ a ( ! ) = 2 R ( ^ R ) : (6.46) Ho w ev er, ^ · ¾ 2 ( ! ) is of little in terest for the follo wing t w o reasons. First, w e also w an t to estimate the phase of ® ( ! ) and hence w e prefer to use (6.45) instead. Second, ev en as an estimate of j ® ( ! ) j 2 , ^ · ¾ 2 ( ! ) in (6.46) can b e sho wn to b e less accurate than the estimate w e can obtain from (6.45) since the latter is obtained b y using the w a v eform structures in (6.9) and (6.15) while the former is not.

PAGE 103

91 6.3.2.2 Determination of ^ a ( ! ) By its v ery de¯nition, ^ a ( ! ) is a v ector in the vicinit y of a ( ! ) suc h that j ® ( ! ) j 2 ^ a ( ! ) ^ a H ( ! ) is a go o d ¯t to ^ R . This leads to the R CB form ulation directly where ^ a ( ! ) is assumed to b elong to a spherical uncertain t y set [38 , 75]: max ¾ 2 ( ! ) ; ^ a ( ! ) ¾ 2 ( ! ) sub ject to ^ R ¡ ¾ 2 ( ! ) ^ a ( ! ) ^ a H ( ! ) ¸ 0 k ^ a ( ! ) ¡ a ( ! ) k 2 · ²; (6.47) where ¾ 2 ( ! ) = j ® ( ! ) j 2 , and ² is a user parameter. A v ector ^ a ( ! ) that is in the range space of ^ S , i.e., ^ a ( ! ) = ^ S ³ for some non-zero ³ , is what w e are after since otherwise w e will get a sp ectral estimate equal to zero, according to (6.44). Hence, for eac h frequency ! , the problem of in terest has the form: max ¾ 2 ( ! ) ; ³ ¾ 2 ( ! ) sub ject to ^ R ¡ ¾ 2 ( ! ) ^ a ( ! ) ^ a H ( ! ) ¸ 0 ^ a ( ! ) = ^ S ³ k ^ a ( ! ) ¡ a ( ! ) k 2 · ²: (6.48) The user parameter ² is used to describ e the uncertain t y of a ( ! ) caused b y the small snapshot n um b er, 2 L < M , in the complex sp ectral estimation problem. The smaller the L , the larger the ² should b e c hosen. Ho w ev er, to exclude the trivial solution of ^ ³ = 0, w e require that ² < k a ( ! ) k 2 = M . Next w e note the follo wing equiv alences for the norm constrain t: k ^ a ( ! ) ¡ a ( ! ) k 2 · ² , ° ° ° ° · ^ S H ^ G H ¸ ( ^ a ( ! ) ¡ a ( ! )) ° ° ° ° 2 · ² , k ^ S H ( ^ a ( ! ) ¡ a ( ! )) k 2 + k ^ G H ( ^ a ( ! ) ¡ a ( ! )) k 2 · ² , k ³ ¡ ¹ ³ k 2 · ² ¡ k ^ G H a ( ! ) k 2 , k ³ ¡ ¹ ³ k 2 · ¹ ²; (6.49) where w e ha v e de¯ned ¹ ³ 4 = ^ S H a ( ! ) (6.50)

PAGE 104

92 and ¹ ² 4 = ² ¡ k ^ G H a ( ! ) k 2 : (6.51) It follo ws from (6.49) that if ¹ ² < 0, whic h o ccurs when a ( ! ) is \far" a w a y from the range space of ^ S , the optimization problem in (6.48) is infeasible: in suc h a case there is no ^ a ( ! ) of the form ^ a ( ! ) = ^ S ³ that satis¯es the constrain t in (6.48), or equiv alen tly , (6.49). Consequen tly , the v ector ^ a ( ! ) cannot b elong to the range space of ^ S in this case and then according to (6.40)-(6.44), w e get ^ ® ( ! ) = 0 : (6.52) Clearly , this indicates that the sp ectral comp onen t at the curren t frequency is negligible. Next, consider the case when ¹ ² ¸ 0, whic h o ccurs when a ( ! ) is \close" to the range space of ^ S and hence of ^ R . In this case, an ^ a ( ! ) b elonging to the range space of ^ R can b e found within the spherical uncertain t y set in (6.49). T o exclude the trivial solution of ³ = 0 (hence ^ a ( ! ) = 0 ), w e assume that k ¹ ³ k 2 > ¹ ²: (6.53) It can b e readily v eri¯ed that the condition in (6.53) is equiv alen t to M = k a ( ! ) k 2 > ² (6.54) (whic h w as assumed b efore). F or an y giv en ^ a ( ! ) in the range of ^ S , the solution ^ ¾ 2 ( ! ) to (6.48) is giv en b y (see App endix 6.6.2): ^ ¾ 2 ( ! ) = 1 ^ a H ( ! ) ^ R y ^ a ( ! ) : (6.55) Making use of the fact that ^ a H ( ! ) ^ R y ^ a ( ! ) = ³ H ^ ª ¡ 1 ³ and of the equiv alences in (6.49), w e can therefore reform ulate (6.48) as the follo wing minimization problem with a quadratic ob jectiv e function and a quadratic inequalit y constrain t: min ³ ³ H ^ ª ¡ 1 ³ sub ject to k ³ ¡ ¹ ³ k 2 · ¹ ²: (6.56)

PAGE 105

93 Because the solution to (6.56) (under (6.53) or (6.54)) will eviden tly o ccur on the b oundary of the constrain t set, w e can reform ulate (6.56) as the follo wing quadratic problem with a quadratic equalit y constrain t: min ³ ³ H ^ ª ¡ 1 ³ sub ject to k ³ ¡ ¹ ³ k 2 = ¹ ²: (6.57) This problem can b e solv ed b y using the L agr ange multiplier metho dolo gy , whic h is based on the function: f = ³ H ^ ª ¡ 1 ³ + ¸ ¡ k ³ ¡ ¹ ³ k 2 ¡ ¹ ² ¢ ; (6.58) where ¸ ¸ 0 is the Lagrange m ultiplier. Di®eren tiation of (6.58) with resp ect to ³ giv es the optimal solution ^ ³ : ^ ª ¡ 1 ^ ³ + ¸ ( ^ ³ ¡ ¹ ³ ) = 0 : (6.59) The ab o v e equation yields ^ ³ = Ã ^ ª ¡ 1 ¸ + I ! ¡ 1 ¹ ³ (6.60) = ¹ ³ ¡ ³ I + ¸ ^ ª ´ ¡ 1 ¹ ³ ; (6.61) where w e ha v e used the matrix in v ersion lemma to obtain the second equalit y . Substituting (6.61) in to the equalit y constrain t of (6.57), the Lagrange m ultiplier ¸ ¸ 0 is obtained as the solution to the constrain t equation: g ( ¸ ) 4 = ° ° ° ° ³ I + ¸ ^ ª ´ ¡ 1 ¹ ³ ° ° ° ° 2 = ¹ ²: (6.62) It can b e sho wn that g ( ¸ ) is a monotonically decreasing function of ¸ ¸ 0 (see, e.g., [38 ]). Hence a unique solution ¸ ¸ 0 exists whic h can b e obtained e±cien tly via, for example, a Newton's metho d. Once ¸ has b een determined, w e use it in (6.61) to get ^ ³ , whic h giv es ^ ¾ 2 ( ! ) = 1 ^ ³ H ^ ª ¡ 1 ^ ³ : (6.63)

PAGE 106

94 Ho w ev er, for the same reasons w e did not use ^ · ¾ 2 ( ! ) in (6.46), w e will not use ^ ¾ 2 ( ! ) as an estimate of j ® ( ! ) j 2 . Instead, w e obtain the rank-de¯cien t R CF h ( ! ) b y substituting ^ a ( ! ) = ^ S ^ ³ , with ^ ³ giv en b y (6.60), in to (6.38) ^ h ( ! ) = ^ S ^ ª ¡ 1 ^ ³ a H ( ! ) ^ S ^ ª ¡ 1 ^ ³ = ^ S ( I ¸ + ^ ª ) ¡ 1 ^ S H a ( ! ) a H ( ! ) ^ S ( I ¸ + ^ ª ) ¡ 1 ^ S H a ( ! ) : (6.64) In App endix 6.6.3, w e sho w that the ^ h ( ! ) deriv ed ab o v e using the ^ R in (6.25) satis¯es J ^ h ¤ ( ! ) = ^ h ( ! ) e ¡ j ( M ¡ 1) ! : (6.65) Then w e ha v e e ¡ j ( N ¡ 1) ! ~ g H ( ! ) ^ h ( ! ) = e ¡ j ( N ¡ 1) ! ¢ 1 L Ã L ¡ 1 X l =0 ~ y l e ¡ j ! l ! H ¢ ^ h ( ! ) = 1 L Ã L ¡ 1 X l =0 ~ y l e ¡ j ! l ! H ¢ J ^ h ¤ ( ! ) ¢ e j ( M ¡ 1) ! e ¡ j ( N ¡ 1) ! = 1 L Ã L ¡ 1 X l =0 J ¹ y ¤ L ¡ l ¡ 1 e ¡ j ! l ! H ¢ J ^ h ¤ ( ! ) ¢ e ¡ j ( N ¡ M ) ! = 1 L Ã J L ¡ 1 X l 0 =0 ¹ y ¤ l 0 e ¡ j ! ( L ¡ l 0 ¡ 1) ! H ¢ J ^ h ¤ ( ! ) ¢ e ¡ j ( N ¡ M ) ! = ( ¹ g ¤ ( ! )) H ¢ e j ( L ¡ 1) ! J ¢ J ^ h ¤ ( ! ) ¢ e ¡ j ( N ¡ M ) ! = ( ¹ g ¤ ( ! )) H ^ h ¤ ( ! ) = ^ h H ( ! ) ¹ g ( ! ) (6.66) Consequen tly , the forw ard-bac kw ard sp ectral estimate ^ ® ( ! ) in (6.45) can b e simpli¯ed as ^ ® ( ! ) = ½ ^ h H ( ! ) ¹ g ( ! ) ; ^ a ( ! ) 2 R ( ^ R ) 0 ; ^ a ( ! ) = 2 R ( ^ R ) ; (6.67) Substituting (6.64) in to (6.67), w e get ^ ® ( ! ) = 8 < : a H ( ! ) ^ S ( I ¸ + ^ ª ) ¡ 1 ^ S H ¹ g ( ! ) a H ( ! ) ^ S ( I ¸ + ^ ª ) ¡ 1 ^ S H a ( ! ) ; ¹ ² ¸ 0 0 ; ¹ ² < 0 : (6.68)

PAGE 107

95 This completes the deriv ation of the rank-de¯cien t R CF sp ectral estimator; a summary of the corresp onding algorithm is giv en b elo w. The selection of ² dep ends on M , L , and the signal-to-noise ratio (SNR). Usually a larger M and a smaller L suggest a larger ² . In practice, w e c ho ose ² exp erimen tally for di®eren t applications based on their di®eren t data size, ¯lter length, n um b er of snapshots, and SNR. W e remark that the ab o v e rank-de¯cien t R CF sp ectral estimator requires O ( K 3 ) °ops, whic h are mainly due to the eigen-decomp osition of the singular sample co v ariance matrix ^ R , while the full rank v ersion needs O ( M 3 ) °ops. 6.4 Numerical Examples W e study the resolution and accuracy p erformance of the rank-de¯cien t R CF complex sp ectral estimator b y using b oth 1-D and 2-D n umerical examples. W e compare the rank-de¯cien t R CF sp ectral estimator with the follo wing sp ectral estimators: the windo w ed FFT (WFFT), Cap on, APES, full-rank norm constrained Cap on ¯lter-bank (NCCF), rank-de¯cien t NCCF, full-rank R CF, and a v ersion of the HDI (see App endix 6.6.5 for brief descriptions of the NCCF and HDI). The v ersion of HDI w e ha v e considered includes b oth norm and subspace constrain ts [3 , 4]. In the ¯rst 1-D example, w e consider estimating the lo cations and complex amplitudes of t w o closely spaced sp ectral lines in the presence of strong in terferences and additiv e zero-mean white Gaussian noise. In the second 1-D example, w e consider a single sp ectral line in the presence of strong in terferences and noise. In the 2-D example, w e in v estigate the usage of complex sp ectral estimators for syn thetic ap erture radar (SAR) imaging. 6.4.1 1-D Complex Sp ectral Estimation W e ¯rst consider t w o closely spaced sp ectral lines (sin usoids) ha ving frequencies 0.09 and 0.1 Hz. F or simplicit y , w e assume that they b oth ha v e unit amplitude

PAGE 108

96 and zero phase. There are 11 strong in terferences that are uniformly spaced b et w een 0.25 and 0.27 Hz in frequency with the frequency spacing b et w een t w o adjacen t interferences b eing 0.002 Hz. The in terferences also ha v e zero phase and are of equal p o w er, whic h is 32 dB stronger than that of the t w o w eak sp ectral lines. The data sequence has 64 samples and is corrupted b y a zero-mean additiv e white Gaussian noise with v ariance ¾ 2 n . F or the t w o sp ectral lines of in terest, w e ha v e the signal-to-noise ratios SNR 1 = SNR 2 = 12 dB, where SNR k = 10 log 10 j ® k j 2 ¾ 2 n (dB) ; (6.69) with ® k b eing the complex amplitude of the k th sin usoid. The true sp ectrum of the signal is giv en in Fig. 6.1(a). W e are in terested in estimating the t w o w eak sp ectral lines. F or b etter visualization, the corresp onding zo omed-in sp ectrum fo cusing on the w eak targets is sho wn in the upp er-left corner of the ¯gure. The mo dulus of the sp ectral estimates obtained b y using WFFT, Cap on, APES, full-rank NCCF, rank-de¯cien t NCCF, full-rank R CF, and rank-de¯cien t R CF are giv en in Figs. 6.1(b)-6.1(h), resp ectiv ely . The comparison with HDI will b e giv en later in a 2-D example. In Fig. 6.1(b), a T a ylor windo w with order 5 and sidelob e lev el -50 dB is applied to the data b efore the zero-padded FFT. Note that the resolution of WFFT is quite p o or. In Figs. 6.1(c) and 6.1(d), resp ectiv ely , Cap on and APES are used, b oth with M = 32. Although the Cap on sp ectrum in Fig. 6.1(c) giv es t w o p eaks close to the desired frequencies, they are not v ery w ell separated. The amplitude estimates of Cap on are also sligh tly biased do wn w ard. The APES sp ectrum is kno wn to giv e excellen t amplitude estimates at the true frequency locations but su®ers from biased frequency estimation [24 ]. As sho wn in Fig. 6.1(d), APES barely resolv es the t w o sp ectral lines. The full-rank NCCF sp ectrum obtained with M = 32 and a norm squared constrain t on the w eigh t v ector corresp onding to ´ = 0 : 3 (see App endix 6.6.5) is sho wn in Fig. 6.1(e); the t w o closely spaced sp ectral lines are hardly separated. The rank-de¯cien t NCCF sp ectrum sho wn in Fig. 6.1(f )

PAGE 109

97 is obtained with M = 56 and a norm squared constrain t of ´ = 2. It has b etter resolution than its full-rank coun terpart in Fig. 6.1(e) but with higher sidelob es, and the amplitude estimates are biased do wn w ard. The full-rank R CF sp ectrum is obtained with M = 32 and ² = 0 : 3, whic h is sho wn in Fig. 6.1(g). Lik e the full-rank NCCF, the full-rank R CF can hardly resolv e the t w o sp ectral lines. Fig. 6.1(h) sho ws the rank-de¯cien t R CF sp ectrum, whic h is obtained b y using M = 56 and ² = 0 : 3. Note that in Fig. 6.1(h), the t w o closely spaced sp ectral lines are w ell resolv ed with no sidelob es. Although APES can pro vide excellen t amplitude estimates at the true frequency lo cations, in man y cases this kno wledge is not a v ailable. When this kno wledge is una v ailable, the frequency estimate for eac h of the t w o sp ectral lines can b e obtained from the cen ter frequency of the half-p o w er (3 dB) in terv al of the corresp onding p eaks in the rank-de¯cien t R CF sp ectrum. Using 100 Mon te Carlo sim ulations (b y v arying the noise realizations), w e obtained the ro ot mean-squared errors (RMSE's) of the frequency estimates of the rank-de¯cien t R CF. F or the ¯rst and second lines of in terest, the RMSE's of the frequency estimates obtained via the rank-de¯cien t R CF are 6 : 3 £ 10 ¡ 4 and 5 : 9 £ 10 ¡ 4 Hz, resp ectiv ely , whic h are quite accurate. The corresp onding RMSE's of the magnitude and phase estimates obtained b y using the rank-de¯cien t R CF and APES at these estimated frequencies are listed in T able I. Note that in this example, the rank-de¯cien t R CF giv es sligh tly more accurate magnitude estimates but w orse phase estimates than APES. T o pro vide further comparisons b et w een the more successful metho ds in the previous example, w e next consider estimating the parameters of a single sp ectral line at 0.09 Hz whic h has unit amplitude and zero phase. The mo dulus of the true complex sp ectrum is sho wn in Fig. 6.2(a). The setup of this exp erimen t is the same as for the previous example except that w e ha v e only one sp ectral line instead of

PAGE 110

98 t w o. Figs. 6.2(b) 6.2(d) sho w the sp ectral estimates obtained with APES, rankde¯cien t NCCF, and rank-de¯cien t R CF, resp ectiv ely . The APES sp ectrum with M = 32 giv es a go o d amplitude estimate and lo w sidelob es, but has a relativ ely wide mainlob e. The rank-de¯cien t NCCF sp ectrum with M = 56 and ´ = 2 sho wn in Fig. 6.2(c) is clearly biased do wn w ard and has high sidelob es. The rank-de¯cien t R CF sp ectrum with M = 56 and ² = 0 : 3 sho wn in Fig. 6.2(d) demonstrates a go o d amplitude estimate, a narro w mainlob e, and no sidelob es. Using 100 Mon te Carlo sim ulations, w e computed the RMSE's of the frequency , magnitude, and phase estimates of the sp ectral line at 0.09 Hz obtained b y using the rank-de¯cien t R CF and APES. The frequency estimates of b oth the rankde¯cien t R CF and APES are obtained b y using the pro cedure for R CF describ ed in the previous example. The RMSE's of the frequency estimates obtained via the rankde¯cien t R CF and APES are listed in T able I I. The RMSE's of the magnitude and phase estimates of the rank-de¯cien t R CF and APES at the frequencies estimated b y the rank-de¯cien t R CF are also listed in T able I I. In this example, the rank-de¯cien t R CF giv es more accurate frequency estimates than APES, but its magnitude and phase estimates are sligh tly w orse. The RMSE's of the magnitude and phase estimates of APES at the frequencies estimated b y APES are 0.028 and 0.057 (radian), resp ectiv ely . These RMSE's are sligh tly w orse than those at the frequencies determined b y the rank-de¯cien t R CF, but they are still sligh tly b etter than those of the rank-de¯cien t R CF estimates. T able 6.1: RMSE's of the mo dulus and the phase (radian) estimates obtained b y the rank-de¯cien t R CF and APES sp ectral estimators in the ¯rst 1-D example. rank-de¯cien t R CF APES mo dulus phase (radian) mo dulus phase (radian) signal 1 0.065 0.393 0.079 0.157 signal 2 0.063 0.415 0.072 0.156

PAGE 111

99 T able 6.2: RMSE's of the frequency , the mo dulus, and the phase (radian) estimates obtained b y the rank-de¯cien t R CF and APES sp ectral estimators in the second 1-D example. rank-de¯cien t R CF APES frequency (Hz) 2 : 72 £ 10 ¡ 4 3 : 15 £ 10 ¡ 4 mo dulus 0.034 0.026 phase (radian) 0.062 0.054 Note that the rank-de¯cien t R CF giv es sidelob e-free sp ectral estimates ( ^ ® ( ! ) = 0 for the case of ¹ ² < 0), whic h also means that w eak signals migh t b e cut o®. (Consider the extreme case where the signal is totally buried in noise.) Clearly , this cuto® e®ect is determined b y SNR and ² . F or the ab o v e 1-D example, the threshold SNR of cutting o® the w eak target is -10 dB for M = 56 and ² = 5. 6.4.2 2-D Complex Sp ectral Estimation No w w e consider using the rank-de¯cien t R CF for SAR imaging. The 2-D high resolution phase history data of a Slicy ob ject at 0 ± azim uth angle w as generated b y XP A TCH [2 ], a high frequency electromagnetic scattering prediction co de for complex 3-D ob jects. A photo of the Slicy ob ject tak en at 45 ± azim uth angle is sho wn in Figure 6.3(a). The original XP A TCH data matrix has a size of N = ¹ N = 288 with a resolution of 0.043 meters in b oth range and cross-range. Fig. 6.3(b) sho ws the mo dulus of the 2-D WFFT image of the original data, where a T a ylor windo w with order 5 and p eak sidelob e lev el -35dB is applied to the data b efore zero-padded FFT. Next, W e consider only a 24 £ 24 cen ter blo c k of the phase history data for SAR image formation, with the purp ose of using Fig. 6.3(b) as a reference for comparison. Since some of the Slicy features, suc h as the sp ectral lines corresp onding to the dihedrals, are not stationary across the cross-range, the in tensit y of the features relativ e to eac h other ma y c hange as the data dimension is reduced from 288 £ 288 to 24 £ 24. Figs. 6.3(c), 6.3(d), 6.3(e), and 6.3(f ) sho w the mo dulus of 2-D FFT,

PAGE 112

100 2-D WFFT (using the same t yp e of windo w as for Fig. 6.3(b)), 2-D Cap on, and 2-D APES sp ectral estimates, resp ectiv ely . Note the high sidelob es and smeared features in the FFT image. The WFFT image demonstrates more smeared features with some of the features not resolv ed. Cap on giv es narro w mainlob es but smaller amplitude estimates than WFFT. APES pro vides un biased sp ectral estimates but has wider mainlob es than Cap on. The mo dulus of the full-rank R CF sp ectral estimate is sho wn in Fig. 6.3(g) whic h is obtained using M = ¹ M = 12 and ² = 2. Fig. 6.3(h) sho ws the mo dulus of the rank-de¯cien t R CF sp ectral estimate obtained b y using M = ¹ M = 16 and ² = 2. Note that the image in Fig. 6.3(h) has no sidelob e problem and all imp ortan t features of the Slicy ob ject are clearly separated. Compared with Fig. 6.3(b), w e note that although the data size w as reduced to 24 £ 24 from 288 £ 288, the rank-de¯cien t R CF pro duces an image similar to the WFFT image using the original high-resolution data. The result for the rank-de¯cien t NCCF is sho wn in Fig. 6.3(i) whic h is obtained using M = ¹ M = 16 and a norm squared constrain t on the w eigh t v ector of ´ = 0 : 2. Compared with Fig. 6.3(h), the features of the rank-de¯cien t NCCF image are not as clear as those of the rank-de¯cien t R CF image and the ¯delit y of the rankde¯cien t NCCF image is w orse. Another rank-de¯cien t sp ectral estimate w e used in this example is the HDI with b oth quadratic and subspace constrain ts [3 ]. The reconstructed image using HDI is sho wn in Fig. 6.3(j) whic h is obtained b y using M = ¹ M = 16 and a norm squared constrain t of % = 0 : 05 (see App endix V). The subspace constrain ts w ere in tro duced in [3 ] to preserv e the bac kground. W e note that the HDI image is more smeared than those obtained using the rank-de¯cien t NCCF and rank-de¯cien t R CF.

PAGE 113

101 6.4.3 Com bining MAPES and R CF In this example, w e apply rank-de¯cien t R CF to the missing-data sp ectral estimation problem where the missing samples are estimated via MAPES. Consider the missing data SAR imaging example in Figure 5.4. Based on the data in terp olated via MAPES-CM, w e apply the rank-de¯cien t robust Cap on ¯lter-bank (R CF) sp ectral estimator with a 20 £ 20 ¯lter-bank and a spherical steering v ector uncertain t y set with unit radius. The resulting image is sho wn in Figure 6.4 (b), whic h exhibits no sidelob e problem and retains all imp ortan t features of Slicy . F or comparison purp ose, the corresp onding 2-D MAPES-CM image is re-plotted in Figure 6.4(a). Compared with Figure 5.4(b), w e note that although the data size w as reduced from 288 £ 288 to 32 £ 32 and from the reduced data matrix 40% of the samples w ere omitted, w e can still obtain an image similar to that obtained b y the WFFT applied to the original high-resolution data. (Note that w e cannot get a similar high resolution image with all the w ell separated features and without sidelob e problem b y simply thresholding Figure 5.4(f ), 5.4(g), or ev en 5.4(c).) 6.5 Conclusions W e ha v e presen ted a new nonparametric FIR ¯ltering based complex sp ectral estimator obtained b y using a rank-de¯cien t robust Cap on ¯lter-bank (R CF) approac h. W e ha v e compared the new sp ectral estimator with other FIR ¯ltering based approac hes including Cap on, APES, NCCF and HDI. W e ha v e found that b y allo wing the sample co v ariance matrix to b e rank-de¯cien t, w e can ac hiev e m uc h higher resolution than the existing sp ectral estimators. W e ha v e sho wn via n umerical examples that for line sp ectral estimation, the frequency estimates obtained b y using the rank-de¯cien t R CF are v ery accurate. The amplitude estimates obtained at these estimated frequencies b y using the rank-de¯cien t R CF are comparable with those obtained b y using APES at the same frequency estimates, whic h are among

PAGE 114

102 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 10 20 30 40 50 Frequency (Hz)Modulus of Complex Amplitude 0 0.1 0.2 0 0.5 1 1.5 ZOOM IN 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 10 20 30 40 50 Frequency (Hz)Modulus of Complex Amplitude 0 0.1 0.2 0 0.5 1 1.5 ZOOM IN (a) (b) 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 10 20 30 40 50 Frequency (Hz)Modulus of Complex Amplitude 0 0.1 0.2 0 0.5 1 1.5 ZOOM IN 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 10 20 30 40 50 Frequency (Hz)Modulus of Complex Amplitude 0 0.1 0.2 0 0.5 1 1.5 ZOOM IN (c) (d) Figure 6.1: Mo dulus of 1-D sp ectral estimates ( N = 64): (a) T rue sp ectrum, (b) WFFT, (c) Cap on with M = 32, (d) APES with M = 32, (e) full-rank NCCF with M = 32 and ´ = 0 : 3, (f ) rank-de¯cien t NCCF with M = 56 and ´ = 2, (g) full-rank R CF with M = 32 and ² = 0 : 3, and (h) rank-de¯cien t R CF with M = 56 and ² = 0 : 3.

PAGE 115

103 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 10 20 30 40 50 Frequency (Hz)Modulus of Complex Amplitude 0 0.1 0.2 0 0.5 1 1.5 ZOOM IN 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 10 20 30 40 50 Frequency (Hz)Modulus of Complex Amplitude 0 0.1 0.2 0 0.5 1 1.5 ZOOM IN (e) (f ) 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 10 20 30 40 50 Frequency (Hz)Modulus of Complex Amplitude 0 0.1 0.2 0 0.5 1 1.5 ZOOM IN 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 10 20 30 40 50 Frequency (Hz)Modulus of Complex Amplitude 0 0.1 0.2 0 0.5 1 1.5 ZOOM IN (g) (h) Figure 6.1: Con tin ued.

PAGE 116

104 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 10 20 30 40 50 Frequency (Hz)Modulus of Complex Amplitude 0 0.1 0.2 0 0.5 1 1.5 ZOOM IN 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 10 20 30 40 50 Frequency (Hz)Modulus of Complex Amplitude 0 0.1 0.2 0 0.5 1 1.5 ZOOM IN (a) (b) 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 10 20 30 40 50 Frequency (Hz)Modulus of Complex Amplitude 0 0.1 0.2 0 0.5 1 1.5 ZOOM IN 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 10 20 30 40 50 Frequency (Hz)Modulus of Complex Amplitude 0 0.1 0.2 0 0.5 1 1.5 ZOOM IN (c) (d) Figure 6.2: Mo dulus of 1-D sp ectral estimates ( N = 64): (a) T rue sp ectrum, (b) APES with M = 32, (c) rank-de¯cien t NCCF with M = 56 and ´ = 2, and (d) rank-de¯cien t R CF with M = 56 and ² = 0 : 3.

PAGE 117

105 -35 -30 -25 -20 -15 -10 -5 0 (a) (b) -35 -30 -25 -20 -15 -10 -5 0 -35 -30 -25 -20 -15 -10 -5 0 (c) (d) Figure 6.3: Mo dulus of the SAR images of the Slicy ob ject from a 24 £ 24 data matrix: (a) Photograph of the ob ject (tak en at 45 ± azim uth angle), (b) 2-D WFFT with 288 £ 288 (not 24 £ 24) data matrix, (c) 2-D FFT, (d) 2-D WFFT, (e) 2-D CAPON with M = ¹ M = 12, (f ) 2-D APES with M = ¹ M = 12, (g) 2-D full-rank R CF with M = ¹ M = 12 and ² = 2, (h) 2-D rank-de¯cien t R CF with M = ¹ M = 16 and ² = 2, (i) 2-D rank-de¯cien t NCCF with M = ¹ M = 16 and ´ = 0 : 2, and (j) 2-D HDI with M = ¹ M = 16 and % = 0 : 05.

PAGE 118

106 -35 -30 -25 -20 -15 -10 -5 0 -35 -30 -25 -20 -15 -10 -5 0 (e) (f ) -35 -30 -25 -20 -15 -10 -5 0 -35 -30 -25 -20 -15 -10 -5 0 (g) (h) -35 -30 -25 -20 -15 -10 -5 0 -35 -30 -25 -20 -15 -10 -5 0 (i) (j) Figure 6.3: Con tin ued.

PAGE 119

107 the b est kno wn. W e ha v e also sho wn via n umerical examples that for 2-D imaging, the rank-de¯cien t R CF sp ectral estimate has the highest ¯delit y as compared with the estimates obtained b y the existing metho ds. The rank-de¯cien t R CF is also an ideal feature extractor due to its prop ert y of placing sp ectral v alues of zero in regions of no in terest. In the pro cess of obtaining the rank-de¯cien t R CF sp ectral estimator, w e ha v e deriv ed a n um b er of results that app ear to b e of indep enden t in terest. In particular, w e ha v e obtained the solutions to b oth the standard form ulation and the co v ariance ¯tting-based form ulation of the Cap on b eamforming problem in the case of a singular sample co v ariance matrix ^ R , and ha v e sho wn that the t w o solutions coincide. Another b y-pro duct of the analysis leading to the R CF sp ectral estimator is a new w a y of motiv ating the robust Cap on approac h (originally in tro duced in [38 ] in a somewhat di®eren t manner): design the ¯lter-bank as in the standard Cap on approac h, y et do not use a but a v ector ^ a in the vicinit y of a that is determined, along with ¾ 2 , so that ¾ 2 ^ a ^ a H pro vides a b etter appro ximation of ^ R (in a certain sense) than ¾ 2 aa H . 6.6 App endixes 6.6.1 2-D Complex Sp ectral Estimation Here w e pro vide a 2-D extension of the forw ard-bac kw ard rank-de¯cien t R CF sp ectral estimator. In the 2-D case, the data mo del in (6.1) can b e written as y n 1 ;n 2 = ® ( ! 1 ; ! 2 ) e j ( ! 1 n 1 + ! 2 n 2 ) + e n 1 ;n 2 ( ! 1 ; ! 2 ) n 1 = 0 ; :::; N 1 ¡ 1; n 2 = 0 ; :::; N 2 ¡ 1; ! 1 ; ! 2 2 [0 ; 2 ¼ ) ; (6.70) where ® ( ! 1 ; ! 2 ) denotes the complex amplitude of a 2-D sin usoidal signal with frequency ( ! 1 ; ! 2 ), and e n 1 ;n 2 ( ! 1 ; ! 2 ) denotes the residual term at frequency ( ! 1 ; ! 2 ). Similar to the 1-D case, w e form the M 1 £ M 2 forw ard-bac kw ard data matrices ¹ Y l 1 ;l 2 = f y ( n 1 ; n 2 ) ; n 1 = l 1 ; :::; l 1 + M 1 + 1 ; n 2 = l 2 ; :::; l 2 + M 2 + 1 g ; ~ Y l 1 ;l 2 = f y ¤ ( n 1 ; n 2 ) ; n 1 = N 1 ¡ l 1 ¡ 1 ; :::; N 1 ¡ l 1 ¡ M 1 ; n 2 = N 2 ¡ l 2 ¡ 1 ; :::; N 2 ¡ l 2 ¡ M 2 g ;

PAGE 120

108 l 1 = 0 ; 1 ; :::; L 1 ¡ 1; l 2 = 0 ; 1 ; :::; L 2 ¡ 1 ; (6.71) where L 1 = N 1 ¡ M 1 + 1 and L 2 = N 2 ¡ M 2 + 1. Let ¹ y l 1 ;l 2 = v ec [ ¹ Y l 1 ;l 2 ] (6.72) ~ y l 1 ;l 2 = v ec [ ~ Y l 1 ;l 2 ] ; (6.73) where v ec [ ¢ ] denotes the op eration of stac king the columns of a matrix on top of eac h other. Let a ( ! 1 ; ! 2 ) = a M 2 ( ! 2 ) ­ a M 1 ( ! 1 ) ; (6.74) where ­ denotes the Kronec k er pro duct, and a M k ( ! k ) = [1 e j ! k ::: e j ( M k ¡ 1) ! k ] T ; k 2 f 1 ; 2 g : (6.75) Then, ¹ y l 1 ;l 2 and ~ y l 1 ;l 2 can b e written as [34 ] ¹ y l 1 ;l 2 = [ ® ( ! 1 ; ! 2 ) a ( ! 1 ; ! 2 )] e j ( ! 1 l 1 + ! 2 l 2 ) + ¹ e l 1 ;l 2 ( ! 1 ; ! 2 ) (6.76) ~ y l 1 ;l 2 = [ ~ ® ( ! 1 ; ! 2 ) a ( ! 1 ; ! 2 )] e j ( ! 1 l 1 + ! 2 l 2 ) + ~ e l 1 ;l 2 ( ! 1 ; ! 2 ) ; (6.77) where ~ ® ( ! 1 ; ! 2 ) = ® ¤ ( ! 1 ; ! 2 ) e ¡ j ( N 1 ¡ 1) ! 1 e ¡ j ( N 2 ¡ 1) ! 2 (6.78) and ¹ e l 1 ;l 2 ( ! 1 ; ! 2 ) and ~ e l 1 ;l 2 ( ! 1 ; ! 2 ) are, resp ectiv ely , formed from f e n 1 ;n 2 ( ! 1 ; ! 2 ) g in the same w a ys as ¹ y l 1 ;l 2 and ~ y l 1 ;l 2 are made from f y n 1 ;n 2 g . The forw ard-bac kw ard sample co v ariance matrix tak es the form of (6.25) where ^ ¹ R and ^ ~ R denotes the sample co v ariance matrices estimated from f ¹ y l 1 ;l 2 g and f ~ y l 1 ;l 2 g , resp ectiv ely , ^ R f = 1 L 1 L 2 L 1 ¡ 1 X l 1 =0 L 2 ¡ 1 X l 2 =0 ¹ y l 1 ;l 2 ¹ y H l 1 ;l 2 (6.79) ^ R b = 1 L 1 L 2 L 1 ¡ 1 X l 1 =0 L 2 ¡ 1 X l 2 =0 ~ y l 1 ;l 2 ~ y H l 1 ;l 2 : (6.80) By making use of the fact that ~ y l 1 ;l 2 = J ¹ y ¤ L 1 ¡ l 1 ¡ 1 ;L 2 ¡ l 2 ¡ 1 ; (6.81)

PAGE 121

109 w e can see that ^ R is also p ersymmetric. Based on the 2-D extension ab o v e, the 2-D rank-de¯cien t R CF sp ectral estimator can b e implemen ted b y follo wing the same steps as summarized in Section 6.3. 6.6.2 Rank-De¯cien t Cap on Beamformer via Co v ariance Fitting In this app endix w e pro v e (6.55) and pro vide some additional insigh ts in to the problem (6.48). Consider the problem (6.48) with a ¯xed ^ a ( ! ), whic h can also b e sho wn to b e a co v ariance matrix ¯tting-based reform ulation of the Cap on b eamformer (see [75 ] and [43 ]): max ¾ 2 ¾ 2 sub ject to ^ R ¡ ¾ 2 ^ a ^ a H ¸ 0 ; (6.82) where ^ a is a giv en v ector, ¾ 2 is the signal p o w er, and ^ R is a singular sample co v ariance matrix. W e solv e the ab o v e optimization problem in (6.82) b y ¯rst considering the case where ^ a b elongs to the range space of ^ R , i.e., ^ a = ^ S ³ ; (6.83) where ³ is a non-zero K £ 1 v ector. F rom the p ositiv e semide¯nite constrain t in (6.82), w e ha v e that: ^ R ¡ ¾ 2 ^ a ^ a H ¸ 0 , ^ S ^ ª ^ S H ¡ ¾ 2 ^ S ³ ³ H ^ S H ¸ 0 , ^ S h ^ ª ¡ ¾ 2 ³ ³ H i ^ S H ¸ 0 , ^ S ^ ª 1 = 2 h I ¡ ¾ 2 ^ ª ¡ 1 = 2 ³ ³ H ^ ª ¡ 1 = 2 i ^ ª 1 = 2 ^ S H ¸ 0 , ¹ S £ I ¡ ¾ 2 ¹ v ¹ v H ¤ ¹ S H ¸ 0 , I ¡ ¾ 2 ¹ v ¹ v H ¸ 0 , 1 ¡ ¾ 2 ¹ v H ¹ v ¸ 0 , ¾ 2 · 1 ¹ v H ¹ v ; (6.84)

PAGE 122

110 where w e ha v e de¯ned ¹ S 4 = ^ S ^ ª 1 = 2 and ¹ v 4 = ^ ª ¡ 1 = 2 ³ . Since ^ a = ^ S ³ , w e ha v e ³ = ^ S H ^ a . Hence ¹ v H ¹ v = ³ H ^ ª ¡ 1 = 2 ^ ª ¡ 1 = 2 ³ = ^ a H ^ S ^ ª ¡ 1 ^ S H ^ a = ^ a H ^ R y ^ a ; (6.85) where ^ R y = ^ S ^ ª ¡ 1 ^ S H is the Mo ore-P enrose pseudo-in v erse of ^ R . Hence the solution to the optimization problem in (6.82) is ^ ¾ 2 = 1 ^ a H ^ R y ^ a : (6.86) whic h pro v es (6.55). Consider no w the case of an ^ a v ector that do es not b elong to the range space of ^ R . Then ^ a can b e written as ^ a = ^ S ³ + ^ G ¯ = £ ^ S ^ G ¤ · ³ ¯ ¸ : (6.87) Let ^ R = £ ^ S ^ G ¤ · ^ ª 0 0 0 ¸ · ^ S H ^ G H ¸ : (6.88) Then w e ha v e that: ^ R ¡ ¾ 2 ^ a ^ a H = £ ^ S ^ G ¤ · ^ ª ¡ ¾ 2 ³ ³ H ¡ ¾ 2 ³ ¯ H ¡ ¾ 2 ¯ ³ H ¡ ¾ 2 ¯ ¯ H ¸ · ^ S H ^ G H ¸ : (6.89) Clearly , if ¯ 6 = 0 , then ^ ¾ 2 = 0 is the only solution to (6.82) for whic h ^ R ¡ ¾ 2 ^ a ^ a H ¸ 0. F or this case, the rank-de¯cien t Cap on b eamformer will giv e an estimated p o w er sp ectrum of zero, i.e., ^ ¾ 2 = 0. Hence the p o w er estimate giv en b y the rank-de¯cien t Cap on b eamformer is: ^ ¾ 2 = ( 1 ^ a H ^ R y ^ a ; ^ a 2 R ( ^ R ) 0 ; ^ a = 2 R ( ^ R ) : (6.90) Note that (6.90) is a scaled v ersion of (6.46), with the scaling factor j ½ j 2 b eing caused b y the constrain t h H ( ! ) ^ a ( ! ) = ½ in (6.29), whic h indicates that the co v ariance

PAGE 123

111 matrix ¯tting form ulation ab o v e is equiv alen t to the standard Cap on b eamforming form ulation ev en in the case of a rank-de¯cien t ^ R . This result is an extension of the one in [43 ] (see also [75 ]), where ^ R w as assumed to b e full-rank. 6.6.3 The Conjugate Symmetry of the F orw ard-Bac kw ard FIR Filter F or the h ( ! ) in (6.64) w e ha v e that: Jh ¤ ( ! ) = J ^ S ¤ ( I ¸ + ^ ª ) ¡ 1 ( ^ S H ) ¤ a ¤ ( ! ) a H ( ! ) ^ S ( I ¸ + ^ ª ) ¡ 1 ^ S H a ( ! ) = J ^ S ¤ ( I ¸ + ^ ª ) ¡ 1 ( ^ S H ) ¤ J ¢ Ja ¤ ( ! ) a H ( ! ) ^ S ( I ¸ + ^ ª ) ¡ 1 ^ S H a ( ! ) = J ^ S ¤ ( I ¸ + ^ ª ) ¡ 1 ( ^ S H ) ¤ J ¢ a ( ! ) ¢ e ¡ j ( M ¡ 1) ! a H ( ! ) ^ S ( I ¸ + ^ ª ) ¡ 1 ^ S H a ( ! ) ; (6.91) where w e ha v e used the facts that a H ( ! ) ^ S ( I ¸ + ^ ª ) ¡ 1 ^ S H a ( ! ) is a real-v alued scalar, and that JJ = I . F rom the p ersymmetry of ^ R , w e ha v e that: J ( ^ S ^ ª ^ S H ) ¤ J = ^ S ^ ª ^ S H , J ^ S ¤ ^ ª ( J ^ S ¤ ) H = ^ S ^ ª ^ S H (6.92) , J ^ S ¤ = ^ S ¢ ; (6.93) where ¢ is a K £ K unitary matrix since ¢ H ¢ = ¢ H ^ S H ^ S¢ = ( ^ S¢ ) H ( ^ S¢ ) = ( J ^ S ¤ ) H ( J ^ S ¤ ) = I : (6.94) Note from (6.92) that w e also ha v e ^ S¢ ^ ª ( J ^ S ¤ ) H = ^ S ^ ª ^ S H ) ¢ ^ ª ( J ^ S ¤ ) H = ^ ª ^ S H ) ¢ ^ ª = ^ ª ^ S H ( J ^ S ¤ ) = ^ ª¢ ; (6.95) whic h implies (along with (6.94)) that ¢ ^ ª¢ H = ^ ª¢¢ H = ^ ª : (6.96)

PAGE 124

112 F rom (6.94) and (6.96), w e get J ( ^ S µ I ¸ + ^ ª ¶ ¡ 1 ^ S H ) ¤ J = ^ S ¢ µ I ¸ + ^ ª ¶ ¡ 1 ¢ H ^ S H = ^ S · ¢ µ I ¸ + ^ ª ¶ ¢ H ¸ ¡ 1 ^ S H = ^ S µ I ¸ + ^ ª ¶ ¡ 1 ^ S H : (6.97) Substituting (6.97) in to (6.91), w e get (6.65). 6.6.4 On the In terpla y Bet w een Sampling Errors and Steering Errors The follo wing discussion is based on [72 , Complemen t 6.5.5]. In general, there is a close relationship b et w een the case of steering v ector errors and small-sample errors, see e.g., [15 ]. More precisely , the sampling estimation errors of the co v ariance matrix can b e view ed as steering v ector errors in a corresp onding theoretical co v ariance matrix, and vice v ersa. F or example, consider a uniform linear arra y (ULA) and assume that the source signals are uncorrelated with one another. In this case, the theoretical co v ariance matrix R of the arra y output is T o eplitz. Assume that the sample co v ariance matrix ^ R is also T o eplitz. According to the Carath ¶ eo dory parameterization of T o eplitz matrices (see [72 , Complemen t 4.9.2]), w e can view ^ R as b eing the theoretical co v ariance matrix asso ciated with a ¯ctitious ULA on whic h uncorrelated signals impinge, but the p o w ers and the direction-of-arriv al's of the latter signals are di®eren t from those of the actual signals. Hence, the small sample estimation errors in ^ R can b e view ed as b eing due to steering v ector errors in a corresp onding theoretical co v ariance matrix. 6.6.5 F orm ulations of NCCF and HDI The norm constrained Cap on ¯lter-bank (NCCF) is giv en b y the solution to the follo wing optimization problem: min h ( ! ) h H ( ! ) ^ R h ( ! ) sub ject to h H ( ! ) a ( ! ) = 1 k h ( ! ) k 2 · ´ ; (6.98)

PAGE 125

113 where the squared norm of the ¯lter is constrained b y ´ . By imp osing this norm constrain t, NCCF can impro v e the robustness of the Cap on ¯lter-bank against the small-sample errors of the co v ariance matrix. Lik e R CF, the resulting NCCF has the form of diagonal loading, where the loading lev el is determined adaptiv ely at di®eren t ! . The subspace constrained and norm constrained HDI ¯lter-bank is giv en b y the solution to the follo wing optimization problem: min h ( ! ) h H ( ! ) ^ R h ( ! ) sub ject to h ( ! ) = a ( ! ) M + º º = ^ S » a H ( ! ) º = 0 k h ( ! ) k 2 · %; (6.99) where the squared norm of the ¯lter is constrained b y % , the subspace constrain t º = ^ S » is used to prev en t h ( ! ) from b ecoming orthogonal to the range space of ^ R , and a H ( ! ) º = 0 is emplo y ed to mak e sure that the solution of º is unique. By adding a subspace constrain t, the HDI ¯lter-bank can preserv e the bac kground (or clutter) regions of the sp ectrum.

PAGE 126

114 (a) (b) Figure 6.4: Mo dulus of the SAR images of the Slicy ob ject obtained from a 32 £ 32 data matrix with missing samples: (a) 2-D MAPES-CM with M 1 = M 2 = 16, (b) 2-D MAPES-CM follo w ed b y 2-D rank-de¯cien t R CF with a 20 £ 20 ¯lter-bank and a unit radius spherical constrain t.

PAGE 127

CHAPTER 7 CONCLUSIONS AND FUTURE RESEAR CH 7.1 Summary of Results In the literature of sp ectral estimation, man y adv anced sp ectral estimation metho ds ha v e b een prop osed, including parametric and nonparametric adaptiv e ¯ltering based approac hes. In general, the nonparametric approac hes are less sensitiv e to data mismo delling than their parametric coun terparts. Moreo v er, the adaptiv e ¯lter-bank based nonparametric sp ectral estimators can pro vide high resolution, lo w sidelob es, and accurate sp ectral estimates while retaining the robust nature of the nonparametric metho ds. This researc h is aimed to in v estigate adaptiv e ¯lter-bank based nonparametric complex sp ectral estimation of data sequences (one-dimensional case) or matrices (t w o-dimensional case) with missing samples o ccurring in arbitrary patterns. The main results of this w ork can b e summarized as follo ws. First, w e consider the one-dimensional (1-D) case. W e dev elop t w o missingdata amplitude and phase estimation (MAPES) algorithms b y using a \maxim um lik eliho o d (ML) ¯tting" criterion. Then w e use the w ell-kno wn exp ectation maximization (EM) metho d to solv e the so-obtained estimation problem iterativ ely . Through n umerical sim ulations, w e demonstrate the excellen t p erformance of the MAPES algorithms for missing-data sp ectral estimation and missing data restoration. Next, w e presen t t w o-dimensional (2-D) extensions of the MAPES-EM algorithms in tro duced in the 1-D case. Then w e dev elop a new MAPES algorithm, referred to as MAPES-CM, b y solving an ML problem iterativ ely via cyclic maximization (CM). W e ha v e compared MAPES-EM with MAPES-CM and ha v e sho wn 115

PAGE 128

116 that MAPES-CM allo ws signi¯can t computational sa vings compared with MAPESEM, whic h is esp ecially desirable for 2-D applications. Numerical examples ha v e b een pro vided to demonstrate the p erformance of the 2-D MAPES algorithms. Finally , b y observing that all the missing-data algorithms dev elop ed previously estimate the missing data samples, w e prop ose to ac hiev e b etter sp ectral estimation p erformance, for example, higher resolution than APES, based on the (complete) data in terp olated via MAPES. It is kno wn that the Cap on sp ectral estimator has b etter resolution (narro w er mainlob es) compared with APES. In this part of the w ork, w e dev elop a rank-de¯cien t robust Cap on ¯lter-bank sp ectral estimator, whic h can ac hiev e m uc h higher resolution than the existing nonparametric sp ectral estimators. 7.2 T opics for F uture Researc h There are sev eral in teresting researc h questions that can b e addressed in the future. In the missing-data sp ectral estimation problem, the follo wing topics are of in terest. ² In Chapters 3 and 5, w e dev elop a group of missing-data algorithms for b oth one-dimensional and t w o-dimensional problems. It is in teresting to in v estigate whether w e can implemen t them more e±cien tly than their presen t forms. ² The same missing-data sp ectral analysis framew ork can b e applied to other sp ectral estimation problems suc h as t w o-dimensional irregularly shap ed data sp ectral analysis. New algorithms could b e dev elop ed b y extending the idea of adaptiv e ¯ltering. ² The missing-data algorithms are dev elop ed based on a \maxim um lik eliho o d (ML) ¯tting" explanation of the APES ¯lter. Since the noise v ectors in our case are o v erlapping, they are not indep enden t of eac h other. Ho w ev er, they are kno wn to ha v e go o d p erformance. More explanations b ey ond \appro ximated ML" migh t exist. The results of rank-de¯cien t R CF sp ectral estimator giv es rise to some relev an t questions. ² It is quite un usual to get completely sidelob e-free sp ectral estimates lik e those sho wn in Chapter 6. Hence w e w ould lik e to kno w theoretically , what is the threshold b elo w that a w eak sp ectral comp onen t of a useful signal migh t b e cut a w a y , as if it w ere wrongly in terpreted as noise. ² The selection of optimal ¯lter length can also b e discussed in the future w ork.

PAGE 129

REFERENCES [1] H.-M. Adorf. In terp olation of irregularly sampled data series { a surv ey . Astr onomic al Data A nalysis Softwar e and Systems IV, ASP Confer enc e Series, V ol. 77 , pages 460{463, 1995. [2] D. J. Andersh, M. Hazlett, S. W. Lee, D. D. Reev es, D. P . Sulliv an, and Y. Ch u. XP A TCH: a high-frequency electromagnetic scattering prediction co de and environmen t for complex three-dimensional ob jects. IEEE A ntennas and Pr op agation Magazine , 36(1):65{69, F ebruary 1994. [3] G. R. Benitz. High-de¯nition v ector imaging. MIT Linc oln L ab or atory JournalSp e cial Issue on Sup err esolution , 10(2):147{170, 1997. [4] G. R. Benitz. High de¯nition v ector imaging for syn thetic ap erture radar. Pr oc e e dings of the 31st Asilomar Confer enc e on Signals, Systems and Computers, P aci¯c Gro v e, CA, pages 1204{1209, No v em b er 1997. [5] R. B. Blac kman and J. W. T uk ey . The me asur ement of p ower sp e ctr a fr om the p oint of view of c ommunic ation engine ering . Do v er Publications, Incorp orated, Do v er, New Y ork, 1959. [6] P . Bro ersen, S. de W aele, and R. Bos. Estimation of autoregressiv e sp ectra with randomly missing data. Pr o c e e dings of the 20th IEEE Instrumentation and Me asur ement T e chnolo gy Confer enc e, V ail, CO, 2:1154{1159, Ma y 2003. [7] T. Bronez. Sp ectral estimation of irregularly sampled m ultidimensional pro cesses b y generalized prolate spheroidal sequences. IEEE T r ansactions on A c oustics, Sp e e ch, and Signal Pr o c essing , 36(12):1862{1873, Decem b er 1988. [8] T. M. Bro wn and J. Christensen-Dalsgaard. A tec hnique for estimating complicated p o w er sp ectra from time series with gaps. The Astr ophysic al Journal , 349:667{674, F ebruary 1990. [9] J. P . Burg. Maxim um en trop y sp ectral analysis. Pr o c e e dings of the 37th Me eting So ciety of Explor ation Ge ophysicists, Oklahoma Cit y , OK, Octob er 1967. [10] J. Cap on. High resolution frequency-w a v en um b er sp ectrum analysis. Pr o c e e dings of the IEEE , 57:1408{1418, August 1969. [11] G. Casella and R. L. Berger. Statistic al Infer enc e . Duxbury , P aci¯c Gro v e, CA, 2001. [12] K. M. Cuomo, J. E. Piou, and J. T. Ma yhan. Ultra wide-band coheren t pro cessing. IEEE T r ansactions on A ntennas and Pr op agation , 47(6):1094{1107, June 1999. [13] A. P . Dempster, N. M. Laird, and D. B. Rubin. Maxim um-lik eliho o d from incomplete data via the EM algorithm. Journal of the R oyal Statistic al So ciety , 39(1):1{38, 1977. 117

PAGE 130

118 [14] G. G. F ahlman and T. J. Ulryc h. A new metho d for estimating the p o w er sp ectrum of gapp ed data. Monthly Notic es, R oyal Astr onomic al So ciety , 199:53{ 65, 1982. [15] D. D. F eldman and L. J. Gri±ths. A pro jection approac h for robust adaptiv e b eamforming. IEEE T r ansactions on Signal Pr o c essing , 42:867{876, April 1994. [16] A. F errari and G. Alengrin. AR sp ectral analysis with random missing observ ations. Pr o c e e dings of the Ninth IEEE Signal Pr o c essing Workshop on Statistic al Signal and A rr ay Pr o c essing , pages 320{323, Septem b er 1998, IEEE, Los Alamitos, CA. [17] I. F o dor and P . Stark. Multitap er sp ectrum estimation for time series with gaps. IEEE T r ansactions on Signal Pr o c essing , 48(12):3472{3483, Decem b er 2000. [18] W. H. F o y . Comparison of metho ds for sp ectral estimation with in terrupted data. IEEE T r ansactions on Signal Pr o c essing , 41(3):1449{1453, Marc h 1993. [19] F. Gini and F. Lom bardini. Multilo ok APES for m ultibaseline SAR in terferometry . IEEE T r ansactions on Signal Pr o c essing , 50(7):1800{1803, July 2002. [20] N. J. Grossbard and E. M. Dew an. Metho ds for estimating the auto correlation and p o w er sp ectral densit y functions when there are man y missing data v alues. The Fifth ASSP Workshop on Sp e ctrum Estimation and Mo deling , pages 30{34, Octob er 1990. [21] J. A. HÄ ogb om. Ap erture syn thesis with a non-regular distribution of in terferometer baselines. Astr onomy and Astr ophysics Supplements , 15:417{426, 1974. [22] J.-C. Hung, B.-S. Chen, W.-S. Hou, and L.-M. Chen. Sp ectral estimation under nature missing data. Pr o c e e dings of ICASSP 01 , 5:3061{3064, Ma y 2001. [23] A. Isaksson. Iden ti¯cation of ARX-mo dels sub ject to missing data. IEEE T r ansactions on A utomatic Contr ol , 38(5):813{819, Ma y 1993. [24] A. Jak obsson and P . Stoica. Com bining Cap on and APES for estimation of sp ectral lines. Cir cuits, Systems, and Signal Pr o c essing , 19(2):159{169, 2000. [25] S. M. Ka y . Mo dern Sp e ctr al Estimation: The ory and Applic ation. Pren tice-Hall, Englew o o d Cli®s, NJ, 1988. [26] S. M. Ka y . F undamentals of Statistic al Signal Pr o c essing: Estimation The ory . Pren tice Hall, Upp er Saddle Riv er, NJ, 1993. [27] J. Kuhn. Reco v ering sp ectral information from unev enly sampled data: t w o mac hine-e±cien t solutions. The Astr onomic al Journal , 87(1):196{202, Jan uary 1982. [28] R. Lacoss. Data adaptiv e sp ectral analysis metho ds. Ge ophysics , 36:661{675, 1971. [29] E. Larsson and J. Li. Sp ectral analysis of p erio dically gapp ed data. IEEE T r ansactions on A er osp ac e and Ele ctr onic Systems , 39(3):1089{1097, July 2003. [30] E. Larsson, G. Liu, P . Stoica, and J. Li. High-resolution SAR imaging with angular div ersit y . IEEE T r ansactions on A er osp ac e and Ele ctr onic Systems , 37(4):1359{1372, Octob er 2001.

PAGE 131

119 [31] E. Larsson and P . Stoica. F ast implemen tation of t w o-dimensional APES and CAPON sp ectral estimators. Multidimensional Systems and Signal Pr o c essing , 13:35{54, Jan uary 2002. [32] E. Larsson, P . Stoica, and J. Li. Amplitude sp ectrum estimation for t w odimensional gapp ed data. IEEE T r ansactions on Signal Pr o c essing , 50(6):1343{ 1354, June 2002. [33] E. G. Larsson, J. Li, and P . Stoica. \High-resolution nonparametric sp ectral analysis: theory and applications." In High-R esolution and R obust Signal Pr oc essing, Y. Hua, A. B. Gershman and Q. Cheng (eds.). Marcel-Dekk er, Inc, New Y ork, NY, 2003. [34] H. Li, J. Li, and P . Stoica. P erformance analysis of forw ard-bac kw ard matc hed¯lterbank sp ectral estimators. IEEE T r ansactions on Signal Pr o c essing , 46:1954{ 1966, July 1998. [35] H. Li, W. Sun, P . Stoica, and J. Li. Tw o-dimensional system iden ti¯cation using amplitude estimation. IEEE Signal Pr o c essing L etters , 9(2):61{63, F ebruary 2002. [36] J. Li and P . Stoica. An adaptiv e ¯ltering approac h to sp ectral estimation and SAR imaging. IEEE T r ansactions on Signal Pr o c essing , 44(6):1469{1484, June 1996. [37] J. Li and P . Stoica. E±cien t mixed-sp ectrum estimation with applications to target feature extraction. IEEE T r ansactions on Signal Pr o c essing , 44:281{295, F ebruary 1996. [38] J. Li, P . Stoica, and Z. W ang. On robust Cap on b eamforming and diagonal loading. IEEE T r ansactions on Signal Pr o c essing , 51(7):1702{1715, July 2003. [39] Z. Liu, H. Li, and J. Li. E±cien t implemen tation of Cap on and APES for sp ectral estimation. IEEE T r ansactions on A er osp ac e and Ele ctr onic Systems , 34(4):1314{1319, Octob er 1998. [40] L. Ljung. System Identi¯c ation: The ory for the User . Pren tice-Hall, Inc., Englew o o d Cli®s, NJ 07632, second edition, 1999. [41] N. Lom b. Least-squares frequency analysis of unequally spaced data. Astr ophysics and Sp ac e Scienc e , pages 447{462, 1976. [42] S. L. Marple, Jr. Digital Sp e ctr al A nalysis with Applic ations . Upp er Saddle Riv er, Pren tice Hall NJ, 1987. [43] T. L. Marzetta. A new in terpretation for Cap on's maxim um lik eliho o d metho d of frequency-w a v en um b er sp ectrum estimation. IEEE T r ansactions on A c oustics, Sp e e ch, and Signal Pr o c essing , 31(2):445{449, April 1983. [44] D. D. Meisel. F ourier transforms of data sampled in unequally spaced segmen ts. The Astr onomic al Journal , 84(1):116{126, Jan uary 1979. [45] S. Mirsaidi and J. Oksman. AR sp ectral estimation with randomly missed observ ations. Pr o c e e dings of the 8th IEEE Signal Pr o c essing Workshop on Statistic al Signal and A rr ay Pr o c essing , pages 52{55, June 1996, IEEE, Los Alamitos, CA.

PAGE 132

120 [46] T. K. Mo on. The exp ectation-maximization algorithm. IEEE Signal Pr o c essing Magazine , pages 47{60, No v em b er 1996. [47] L. M. No v ak, G. J. Owirk a, and A. L. W ea v er. Automatic target recognition using enhanced resolution SAR data. IEEE T r ansactions on A er osp ac e and Ele ctr onic Systems , 35(1):157{175, Jan uary 1999. [48] M. R. P alsetia and J. Li. Using APES for in terferometric SAR imaging. IEEE T r ansactions on Imaging Pr o c essing , 7(9):1340{1353, Septem b er 1998. [49] V. F. Pisarenk o. The retriev al of harmonics from a co v ariance function. Ge ophysic al Journal of the R oyal Astr onomic al So ciety , 33(3):347{366, 1973. [50] B. P orat and B. F riedlander. ARMA sp ectral estimation of time series with missing observ ations. IEEE T r ansactions on Information The ory , 30(4):601{ 602, July 1986. [51] W. Press and G. Rybic ki. F ast algorithm for sp ectral analysis of unev enly sampled data. The Astr ophysic al Journal , pages 277{280, Marc h 1989. [52] R. A. Redner and H. F. W alk er. Mixture densities, maxim um lik eliho o d and the EM algorithm. SIAM R eview , 26(2):195{239, April 1984. [53] D. H. Rob erts, J. Leh¶ ar, and J. W. Dreher. Time series analysis with CLEAN I: Deriv ation of a sp ectrum. The Astr onomic al Journal , 93(4):968{989, April 1987. [54] Y. Rosen and B. P orat. Optimal ARMA parameter estimation based on the sample co v ariances for data with missing observ ations. IEEE T r ansactions on Information The ory , 35(2):342{349, Marc h 1989. [55] Y. Rosen and B. P orat. The second-order momen ts of the sample co v ariances for time series with missing observ ations. IEEE T r ansactions on Information The ory , 35(2):334{341, Marc h 1989. [56] R. Ro y , A. P aulra j, and T. Kailath. ESPRIT { A subspace rotation approac h to estimation of parameters of sin usoids in noise. IEEE T r ansactions on A c oustics, Sp e e ch, and Signal Pr o c essing , ASSP-34(5):1340{1342, Octob er 1986. [57] N. Rozario and A. P ap oulis. Sp ectral estimation from nonconsecutiv e data. IEEE T r ansactions on Information The ory , 33(6):889{894, No v em b er 1987. [58] D. J. Russell and R. D. P almer. Application of APES to adaptiv e arra ys on the CDMA rev erse c hannel. IEEE T r ansactions on V ehicular T e chnolo gy , 53(1):3{ 17, Jan uary 2004. [59] G. B. Rybic ki and W. H. Press. In terp olation, realization, and reconstruction of noisy , irregularly sampled data. The Astr ophysic al Journal , 398:169{176, Octob er 1992. [60] J. Salzman, D. Ak amine, and R. Lefevre. Optimal w a v eforms and pro cessing for sparse frequency UWB op eration. Pr o c e e dings of the 2001 IEEE R adar Conferenc e, A tlan ta, GA, pages 105{110, Ma y 2001. [61] J. Salzman, D. Ak amine, R. Lefevre, and J. Kirk, Jr. In terrupted syn thetic ap erture radar (SAR). Pr o c e e dings of the 2001 IEEE R adar Confer enc e, A tlan ta, GA, pages 117{122, Ma y 2001.

PAGE 133

121 [62] J. D. Scargle. Studies in astronomical time series analysis I I: Statistical asp ects of sp ectral analysis of unev enly spaced data. The Astr ophysic al Journal , 263:835{ 853, Decem b er 1982. [63] J. D. Scargle. Studies in astronomical time series analysis I I I: Fourier transforms, auto correlation functions, and cross-correlation functions of unev enly spaced data. The Astr ophysic al Journal , 343:874{887, August 1989. [64] R. Sc hmidt. A Signal Subsp ac e Appr o ach to Multiple Emitter L o c ation and Sp e ctr al Estimation . Ph.D. Dissertation, Stanford Univ ersit y , CA, No v em b er, 1981. [65] U. J. Sc h w arz. Mathematical-statistical description of the iterativ e b eam remo ving tec hnique (Metho d CLEAN). Astr onomy and Astr ophysics , 65:345{356, 1978. [66] P . Stoica, A. Jak obsson, and J. Li. Sin usoid parameter estimation in the colored noise case: asymptotic Cram ¶ er-Rao b ound, maxim um lik eliho o d and nonlinear least-squares. IEEE T r ansactions on Signal Pr o c essing , 45(8):2048{2059, August 1997. [67] P . Stoica, A. Jak obsson, and J. Li. Matc hed-¯lter bank in terpretation of some sp ectral estimators. Signal Pr o c essing , 66(1):45{59, April 1998. [68] P . Stoica, E. G. Larsson, and J. Li. Adaptiv e ¯lterbank approac h to restoration and sp ectral analysis of gapp ed data. The Astr onomic al Journal , 120:2163{2173, Octob er 2000. [69] P . Stoica, H. Li, and J. Li. A new deriv ation of the APES ¯lter. IEEE Signal Pr o c essing L etters , 6(8):205{206, August 1999. [70] P . Stoica, H. Li, and J. Li. Amplitude estimation for sin usoidal signals with application to system iden ti¯cation. to app e ar in IEEE T r ansactions on Signal Pr o c essing , 48(2), F ebruary 2000. [71] P . Stoica and R. L. Moses. Intr o duction to Sp e ctr al A nalysis. Pren tice-Hall, Englew o o d Cli®s, NJ, 1997. [72] P . Stoica and R. L. Moses. Sp e ctr al A nalysis of Signals . Pren tice-Hall, Englew o o d Cli®s, NJ, 2004, forthcoming. [73] P . Stoica and A. Nehorai. Statistical analysis of t w o nonlinear least-squares estimators of sine-w a v e parameters in the colored-noise case. Cir cuits, Systems, and Signal Pr o c essing , 8(1):3{15, 1989. [74] P . Stoica and Y. Selen. Cyclic minimizers, ma jorization tec hniques, and the exp ectation-maximization algorithm: a refresher. IEEE Signal Pr o c essing Magazine , 21(1):112{114, Jan uary 2004. [75] P . Stoica, Z. W ang, and J. Li. Robust Cap on b eamforming. IEEE Signal Pr oc essing L etters , 10(6):172{175, June 2003. [76] P . Sw an. Discrete Fourier transforms of non uniformly spaced data. The Astr onomic al Journal , 87(11):1608{1615, No v em b er 1982. [77] H. L. V an T rees. Optimum A rr ay Pr o c essing, Part IV of Dete ction, Estimation, and Mo dulation The ory . John Wiley and Sons, Inc., New Y ork, NY, 2002.

PAGE 134

122 [78] Y. W ang, J. Li, and P . Stoica. Rank-de¯cien t robust Cap on ¯lter-bank approac h to complex sp ectral estimation. IEEE T r ansactions on Signal Pr o c essing , In press. [79] Y. W ang, P . Stoica, and J. Li. Tw o-dimensional nonparametric sp ectral analysis in the missing data case. IEEE T r ansactions on A er osp ac e and Ele ctr onics Systems , submitted 2004. [80] Y. W ang, P . Stoica, J. Li, and T. Marzetta. Nonparametric sp ectral analysis with missing data via the EM algorithm. Digital Signal Pr o c essing , In press. [81] P . D. W elc h. The use of fast Fourier transform for estimation of p o w er sp ectra: A metho d based on time a v eraging o v er short, mo di¯ed p erio dograms. IEEE T r ansactions on A udio and Ele ctr o ac oustics , 15:70{76, June 1967. [82] R. W u, Z.-S. Liu, and J. Li. Time-v arying complex sp ectral analysis via recursiv e APES. IEE Pr o c e e dings|R adar, Sonar, and Navigation , 145(6):354{360, Decem b er 1998. [83] P . Zhou and A. D. P oularik as. AR sp ectral estimation of segmen ted time signals. Twenty-Se c ond Southe astern Symp osium on System The ory , pages 536{539, Marc h 1990, Co ok eville, TN.

PAGE 135

BIOGRAPHICAL SKETCH Y an w ei W ang w as b orn in Beijing, China, in 1973. He receiv ed the B.Sc degree from Beijing Univ ersit y of T ec hnology , China, in 1997 and the M.Sc degree from the Univ ersit y of Florida, Gainesville, in 2001, b oth in electrical engineering. Since Jan uary 2000, he has b een a researc h assistan t with the Departmen t of Electrical and Computer Engineering, Univ ersit y of Florida, where he is curren tly pursuing the Ph.D. degree in electrical engineering. His curren t researc h in terests include sp ectral estimation, tomographic imaging, and radar/arra y signal pro cessing. He is a studen t mem b er of IEEE. 123