• TABLE OF CONTENTS
HIDE
 Title Page
 Acknowledgement
 Table of Contents
 Abstract
 Introduction
 Robust capon beamforming (RCB)
 Application of RCB to aeroacou...
 Application of RCB to ultrason...
 Conclusions and future work
 Appendix
 Reference
 Biographical sketch
 Copyright






Title: Robust adaptive array processing
CITATION THUMBNAILS PAGE IMAGE ZOOMABLE
Full Citation
STANDARD VIEW MARC VIEW
Permanent Link: http://ufdc.ufl.edu/UF00082338/00001
 Material Information
Title: Robust adaptive array processing theory and applications
Physical Description: vii, 108 leaves : ill. ; 29 cm.
Language: English
Creator: Wang, Zhisong
Publication Date: 2005
 Subjects
Subject: Electrical and Computer Engineering thesis, Ph. D   ( lcsh )
Dissertations, Academic -- Electrical and Computer Engineering -- UF   ( lcsh )
Genre: bibliography   ( marcgt )
theses   ( marcgt )
non-fiction   ( marcgt )
 Notes
Thesis: Thesis (Ph. D.)--University of Florida, 2005.
Bibliography: Includes bibliographical references.
Statement of Responsibility: by Zhisong Wang.
General Note: Printout.
General Note: Vita.
 Record Information
Bibliographic ID: UF00082338
Volume ID: VID00001
Source Institution: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: aleph - 003320266

Table of Contents
    Title Page
        Page i
        Page ii
    Acknowledgement
        Page iii
    Table of Contents
        Page iv
        Page v
    Abstract
        Page vi
        Page vii
    Introduction
        Page 1
        Page 2
        Page 3
        Page 4
        Page 5
        Page 6
        Page 7
        Page 8
        Page 9
    Robust capon beamforming (RCB)
        Page 10
        Page 11
        Page 12
        Page 13
        Page 14
        Page 15
        Page 16
        Page 17
        Page 18
        Page 19
        Page 20
        Page 21
        Page 22
        Page 23
        Page 24
        Page 25
        Page 26
        Page 27
        Page 28
        Page 29
        Page 30
        Page 31
        Page 32
        Page 33
        Page 34
        Page 35
        Page 36
        Page 37
        Page 38
    Application of RCB to aeroacoustics
        Page 39
        Page 40
        Page 41
        Page 42
        Page 43
        Page 44
        Page 45
        Page 46
        Page 47
        Page 48
        Page 49
        Page 50
        Page 51
        Page 52
        Page 53
        Page 54
        Page 55
        Page 56
    Application of RCB to ultrasonics
        Page 57
        Page 58
        Page 59
        Page 60
        Page 61
        Page 62
        Page 63
        Page 64
        Page 65
        Page 66
        Page 67
        Page 68
        Page 69
        Page 70
        Page 71
        Page 72
        Page 73
        Page 74
        Page 75
        Page 76
        Page 77
        Page 78
        Page 79
        Page 80
        Page 81
        Page 82
        Page 83
        Page 84
        Page 85
        Page 86
        Page 87
        Page 88
        Page 89
    Conclusions and future work
        Page 90
        Page 91
        Page 92
    Appendix
        Page 93
        Page 94
        Page 95
        Page 96
        Page 97
        Page 98
        Page 99
    Reference
        Page 100
        Page 101
        Page 102
        Page 103
        Page 104
        Page 105
        Page 106
        Page 107
    Biographical sketch
        Page 108
        Page 109
        Page 110
    Copyright
        Copyright 1
        Copyright 2
Full Text











ROBUST ADAPTIVE ARRAY PROCESSING: THEORY AND APPLICATIONS


By

ZHISONG WANG















A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY

UNIVERSITY OF FLORIDA


2005



























To my parents and my wife.














ACKNOWLEDGMENTS

I would like to sincerely thank my advisor (Prof. Jian Li) for her support,

encouragement, inspiration, and patience in this research. I am greatly indebted

to her for introducing me to the area of robust adaptive array signal processing,

and for her guidance throughout the development of this dissertation. My special

appreciation is due to Prof. Toshikazu Nishida, Prof. Mark Sheplak, and Prof. Louis

N. Cattafesta III for serving on my supervisory committee, and for their valuable

guidance and constructive comments. I am deeply grateful to Prof. Petre Stoica for

his advice, comments, and suggestions.

I am deeply thankful to my parents (Lin Wang and Yuyue Li), my sisters, and

my brothers-in-law for their endless love, constant encouragement, and support. I

sincerely thank my wife (Yu Chen) for her love, inspiration, patience, and care. I am

especially thankful to my father who was unable to see the completion of my study.

I gratefully acknowledge Dr. Jianhua Liu, Dr. Guoqing Liu, Dr. Xi Li, and

Dr. Renbiao Wu for their help during this work. I wish to thank other fellow graduate

students in the SAL group with whom I had the great pleasure of interacting. Finally,

I would like to thank all of the people who helped me during my Ph.D. study.















TABLE OF CONTENTS

page

ACKNOWLEDGMENTS ............................. iii

ABSTRACT ................... ............... vi

CHAPTER

1 INTRODUCTION .................... ........... 1

1.1 Background on Array Signal Processing .. ............... 1
1.1.1 Data Model ................... ......... 1
1.1.2 Uniform Linear Array .............. .... .... 4
1.1.3 Classification .. .... .... ... .. .. .. .. .. ... 5
1.1.4 Beamforming ................... ........ 6
1.2 Motivation for Robust Adaptive Beamforming ............. 7
1.3 Scope of the Work ................... ......... 9

2 ROBUST CAPON BEAMFORMING (RCB) ................ 10

2.1 Problem Formulation .......................... 10
2.2 Standard Capon Beamforming . . . . . . . .. . ... 13
2.2.1 Spatial Filtering SCB ....................... 13
2.2.2 Covariance Fitting SCB ..................... 14
2.3 Robust Capon Beamforming .. .......... .......... 15
2.3.1 Non-Degenerate Ellipsoidal Uncertainty Set . . . ... 15
2.3.2 Flat Ellipsoidal Uncertainty Set . . . . . . ..... 22
2.4 Diagonal Loading Interpretation of RCB . . . . . . ... 25
2.4.1 Nondegenerate Ellipsoidal Uncertainty Set . . . . ... 26
2.4.2 Flat Ellipsoidal Uncertainty Set . . . . . . ..... 26
2.5 Numerical Examples ................... ........ 27

3 APPLICATION OF RCB TO AEROACOUSTICS . . . . .... 39

3.1 Introduction . . . . . . . . . . . . . . .. 39
3.2 Data Model and Problem Formulation of Acoustic Imaging . . 41
3.3 Constant-Powerwidth RCB . . . . . . . ..... ....... 44
3.4 Constant-Beamwidth RCB . . . . . . . ..... ...... 45
3.5 Numerical Examples ................... ........ 46

4 APPLICATION OF RCB TO ULTRASONICS . . . . . .... 57

4.1 Introduction . . . . . . .. . . . . . . .. . 57
4.2 Problem Formulation ............... . .......... 60
4.3 Time-Delay Based RCB ......................... 64









4.4 Time-Reversal Based RCB ................... ..... 68
4.5 Simulated and Experimental Examples . . . . . . .... 72

5 CONCLUSIONS AND FUTURE WORK . . . . . . . .... 90

5.1 Conclusions .... .... .... .. .. .. ... .. .. ... .. 90
5.2 Future W ork ................... ............. 91

APPENDIX

A LINKING RCB TO VOROBYOV AND COLLEAGUES' APPROACH .93

B CALCULATING THE STEERING VECTOR ................ 96

C LINK BETWEEN RCB AND LORENZ AND BOYD'S APPROACH ... 97

REFERENCES .................... ............... 100

BIOGRAPHICAL SKETCH ................... ........ 108














Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy

ROBUST ADAPTIVE ARRAY PROCESSING: THEORY AND APPLICATIONS

By

Zhisong Wang

May 2005

Chair: Jian Li
Major Department: Electrical and Computer Engineering

Adaptive beamformers have better resolution and much better interference

rejection capability than the data-independent beamformers, provided that the array

steering vector corresponding to the signal of interest is accurately known. However,

the adaptive beamforming techniques rely significantly on the assumptions of the

incident signals, the sensor array, and the surroundings; and are very sensitive to the

inevitable mismatch between the assumptions and the actual characteristics. Conse-

quently, adaptive beamformers often suffer from severe performance degradations in

practical applications.

We developed a novel robust Capon beamformer (RCB) that can maintain

both the robustness of the data-independent beamformers and the adaptivity of the

standard Capon beamformer (SCB) (a type of the adaptive beamformer). We showed

that RCB belongs to the class of diagonal loading approaches. However, the amount

of diagonal loading can be precisely calculated based on the uncertainty set of the

array steering vector. We showed how the proposed robust Capon beamformer can

be efficiently computed at a cost comparable to that of SCB. We also provided deep

insights into the relationships among the recent three robust adaptive beamformers.








Although the robust adaptive beamformers obtained from the two different formula-

tions appear to be rather different, we found that they are equivalent in terms of the

signal-to-interference-plus-noise ratio (SINR). However, RCB also has some unique

features. By comparing these three beamformers for degenerate vs non-degenerate el-

lipsoidal uncertainty sets of the steering vector, we showed the impressive advantages

of RCB over the other two methods: simpler implementation, lower computational

complexity, and more accurate power estimation. In addition, we applied RCB to

various applications, and extended RCB to accommodate their specific goals and

requirements. For the application of aero-acoustic imaging, we devised a constant-

beamwidth RCB and a constant-powerwidth RCB for consistent acoustic imaging.

This means that for an acoustic wideband monopole source with a flat spectrum,

the acoustic image for each frequency bin stays approximately the same. For ul-

trasound imaging, we developed a time-delay based RCB and a time-reversal based

RCB. We demonstrated the excellent performance of RCB and its various extensions

by simulated and experimental examples.















CHAPTER 1
INTRODUCTION

Array signal processing [1]-[8] deals with the problem of extracting information

(the number of sources, the signal waveforms, their spatial locations, etc.) from the

measurements taken from an array of spatially distributed sensors in the propagating

wavefield. The wavefield can be electromagnetic, acoustic, or seismic. Sensor ar-

rays have been extensively used in many diverse applications including radar, sonar,

acoustics, communications, speech processing, seismology, astronomy, and medical

imaging.

The main fields of array signal processing include source localization, source

waveform estimation, source intensity estimation, and channel estimation. The objec-

tives for different applications are not the same. For example, in radar, the locations

of the sources are the most important. In communications, the information-bearing

waveforms are most critical. In acoustic measurements, people are concerned with

the sound pressure levels or intensity levels of the sources. In seismic applications,

the interest lies in the seismic structure, which can be determined by examining the

the space-time propagation effects between the sources and the array. In ultrasound

imaging, the reconstructed tissue characteristics and image contrast are crucial for

diagnosing illnesses.

1.1 Background on Array Signal Processing
1.1.1 Data Model

Because of the different configurations, characteristics, and goals in the dis-

tinct applications, the associated array signal processing problems can be very differ-

ent. We use a narrowband far-field data model to show some concepts of array signal









processing. Wideband signals in the near-field are discussed later. A signal is nar-

rowband if the propagation time across the array is small compared to the inverse of

its bandwidth. In [9], the following three conditions are given to determine whether

a source is in the far-field.
2D2
r > (1.1)

r > D (1.2)

r > A (1.3)

where r is the range from the source, and D is the aperture of the array.

Consider an array comprising M sensors. Assume that there are J narrow-

band sources in the far-field. We assume that the data and weights are complex-

valued, since in many applications a quadrature receiver is used to generate in-phase

and quadrature signals, which are combined to yield the complex baseband signals.

Assume that the output of each sensor contains the superposition of the J planar

waveforms and additive sensor noise. Then the output of the ith sensor is given by
J
Xm(t) = Z am(0j)sj(t) + nm(t) (1.4)
j=1

where sj(t) denotes the jth emitter signal, Oj is the direction of arrival (DOA) of

sj(t), am(Oj) is a complex scaler representing the propagating delay of sj(t) and the
gain and phase adjusted by the mth sensor, and nm(t) is the additive noise. In matrix

notation, we have

x(t) = [a(01) ... a(0j)] s(t) + n(t) = A(0)s(t) + n(t) (1.5)

where x(t) = [xi(t) .. XM(t)]T is the data vector at the time t, s(t) = [sl(t) ... s(t)]T

is the signal vector, n(t) = [ni(t) ... nM(t)]T is the noise vector, 0 is a J-dimensional

parameter vector corresponding to the true DOAs, and a(0O) = [al(0j) .. aM(0j)]T
is the array steering vector representing the array response to a unit wavefront from









the direction 9j. The collection of the array steering vectors over the parameter space
of interest, {a(0j)91j e}, is often called the array manifold.
The array steering vector a(0j) can be written as

a(03) = [Hi(wc)e-iw ... HM(wc)e-icrM]T (1.6)

where Hm(wc) is a complex scalar representing the gain and phase at the mth sensor,

Tm denotes the time needed for the wave to travel from a reference point to the ith
sensor, wc is the center or carrier frequency of the signal. If the sensors are identical
and omnidirectional over the DOA range of interest, and the first sensor of the array
is chosen as the reference point, we can obtain the following simplified form for the

steering vector:

a(Oj) = [1 e-iwcr2 .. e-iwTM]T (1.7)

The theoretical covariance matrix of the array output vector is defined as

R = E{x(t)x*(t)} (1.8)

where E{.} is the expectation operator. In practice, the theoretical covariance matrix
R is not known. Hence a sample covariance matrix needs to be estimated based on
the available data samples. Assume that the array output is sampled at N time
instants, and these snapshots are collected in the columns of an M x N data matrix

XN,
XN = [x(1) ... x(N)] = A(O)SN + NN (1.9)

Then the sample covariance matrix can be obtained by

R= -NXX* (1.10)

which is the maximum likelihood estimate of the theoretical covariance matrix. To
make R nonsingular, at least M snapshots are required. The more snapshots, the
better the sample covariance matrix approximates the theoretical covariance matrix.









The cost, however, is longer data acquisition time and more computation. Given only

a finite number of snapshots, undesirable phenomena such as beampattern distortion

and SNR loss may occur.

1.1.2 Uniform Linear Array

If an array consists of several identical sensors uniformly spaced on a line,

it is called a uniform linear array (ULA). Let d denote the distance between two

neighboring sensors, and let 0 denote the DOA of the signal counterclockwise from

the broadside of the array. If we choose the first sensor of the ULA as the reference

point, we have

S= (m -1)(1.11)
c
where c is the wave propagating velocity. Let A be the signal wavelength. Then we

have
c W,
A = f = (1.12)
fc 27
Define

f = fdsin (1.13)
c
dsin9
ws = 27rfs = wc- (1.14)
c
where w8 is referred to as spatial frequency. Inserting (1.14) into (1.11), and then into

(1.7) gives

a(w,) = [1 e-iW ... e-i(M-1)w]T (1.15)

Note that the steering vector in (1.15) is directly analogous to the tapped delay line

vector used for an FIR filter. In fact, a beamformer based on the data collected at

an array of sensors is a spatial filter analogous to a temporal FIR filter based on the

data collected at one sensor. The vector a(ws) is uniquely defined if and only if the

following condition is satisfied:


ws| < 7r


(1.16)






5

which is equivalent to

dl sino 0 A/2 (1.17)

which is satisfied if

d < A/2 (1.18)

For a ULA, if the intersensor spacing d is larger than A/2, spatial aliasing (grating
lobes in the beampattern) will occur. On the other hand, the spatial resolution of a

beamformer is inversely proportional to the array aperture. Hence, it is conventional

to choose

d = A/2 (1.19)

For the sake of simplicity, we consider a ULA with half-wavelength intersensor spacing

in the numerical examples in Chapter 3. The algorithms, however, can be applied to

arbitrary arrays.

1.1.3 Classification

The array signal processing methods are of two kinds: nonparametric and

parametric. The nonparametric method includes delay-and-sum beamforming [4],

[5], standard Capon beamforming [10]-[12], robust adaptive beamforming (e.g., the
robust Capon beamforming approaches [13]-[15]). The parametric method includes

MUSIC [16]-[23], ESPRIT [24]-[28], maximum likelihood (ML) [29]-[32], and weighted
subspace fitting (WSF) [33], [34]. The delay-and-sum beamforming, standard Capon

beamforming, robust Capon beamforming, and MUSIC are discussed later.

By definition, a nonparametric method makes no assumption of the covariance

structure of the data. It only uses the array covariance matrix and the array steering

vector. In contrast, a parametric method requires knowledge of the model in (1.5).
Furthermore, the noise n(t) is assumed to be spatially white, with its components

having identical variances. Then we obtain


E{n(t)n"(t)} = oI


(1.20)









where ua is the noise power (variance). In addition, the signal covariance matrix

S = E{s(t)s*(t)} (1.21)

is often assumed to be nonsingular, which means that the signals may be partially

correlated but not coherent. Then the theoretical covariance matrix of the array

output vector can be written as

R = E{x(t)x*(t)} = A(0)SA*(0) + a2I (1.22)

If the narrowband signals are uncorrelated, S is a diagonal matrix and
J
R= caa(0,)a*(0j) + 4I (1.23)
j=1

where ({a3j}Uj) are the powers of the J signals impinging on the array. In the

presence of coherent signal and interference, S is a singular matrix.

1.1.4 Beamforming

A beamformer is a processor that applies a weight vector on the output of an

array of sensors to estimate the signal of interest in the presence of noise and inter-

ferences. The goal of beamforming is to perform spatial filtering, passing the signal

of interest from a given direction (or location), and at the same time, suppressing the

background noise and directional interference [35]. Beamforming is an essential task

in array signal processing with ubiquitous applications. Beamforming can be used

for either transmitting or receiving signals. We considered beamforming only for the

purpose of reception.

As its name implies, a data-independent beamformer selects its weight vector

regardless of the incoming data. On the other hand, an adaptive (data-dependent)

beamformer takes advantage of the incoming data and adjusts its weight vector as

a function of the data. Standard data-independent beamformers include the delay-

and-sum approach and methods based on various data-independent weight vectors for









sidelobe control [5], [8]. Adaptive beamformers include the standard Capon beam-

former, linearly constrained minimum variance beamformer, and their extensions.

Adaptive beamformers use the block mode or online mode to update the weight

vectors. The beamformers that exploit a known theoretical covariance matrix are

referred to as optimum beamformers.

1.2 Motivation for Robust Adaptive Beamforming

The data-dependent Capon beamformer adaptively selects the weight vector

to minimize the array output power, subject to the linear constraint that the signal-of-

interest (SOI) does not suffer from any distortion [10], [11]. The Capon beamformer

has better resolution and much better interference rejection capability than the data-

independent beamformer, provided that the array steering vector corresponding to

the SOI is accurately known.

However, in practice, the knowledge of the SOI steering vector is often im-

precise, because of differences between the assumed signal arrival angle and the true

arrival angle; or between the assumed array response and the true array response (ar-

ray calibration errors). Whenever this happens, the Capon beamformer may suppress

the SOI as an interference, which results in significantly underestimated SOI power

and drastically reduced array output signal-to-interference-plus-noise ratio (SINR).

Then the performance of the Capon beamformer may become worse than that of the

standard data-independent beamformers [36], [37].

The same happens when the number of snapshots is relatively small (i.e., about

the same as, or smaller than, the number of sensors). In fact, a close relationship

exists between the cases of steering vector errors and small-sample errors [38] in the

sense that the difference between the sample covariance matrix R (estimated from a

finite number of snapshots) and the corresponding theoretical covariance matrix R

can be viewed as due to steering vector errors.









Many approaches have been proposed during the past 3 decades to improve the

robustness of the Capon beamformer. The literature on robust adaptive beamforming

is extensive [8], [13], [15], [35], [39]-[45].

To account for array steering vector errors, additional linear constraints (in-

cluding point and derivative constraints) can be imposed to improve the robustness

of the Capon beamformer [46]-[49]. However, these constraints are not explicitly re-

lated to the uncertainty of the array steering vector. Moreover, for every additional

linear constraint imposed, the beamformer loses one degree of freedom (DOF) for

interference suppression. These constraints belong to the class of covariance matrix

tapering approaches [50].

Subspace based adaptive beamforming methods [38], [51] require knowledge

of the noise covariance matrix. Hence they are sensitive to imprecise knowledge of

the noise covariance matrix and also to the array steering vector errors. Making

these methods robust against array steering vector errors will not cure them of being

sensitive to imprecise knowledge of the noise covariance matrix.

Diagonal loading (including its extended versions) has been a popular and

widely used approach to improve the robustness of the Capon beamformer [2], [51]-
[69]. One representative of the diagonal loading based approaches is the norm con-
strained Capon beamformer (NCCB), which uses a norm constraint on the weight

vector to improve robustness against array steering vector errors, and control white

noise gain [2], [52]-[55]. However, for NCCB and most other diagonal loading meth-

ods, it is not clear how to choose the diagonal loading level based on information

about the uncertainty of the array steering vector.

Recently some methods with a clear theoretical background [13]-[15], [44],

[45], [61], [62] (unlike early methods) make explicit use of an uncertainty set of the
array steering vector. In [61], a polyhedron is used to describe the uncertainty set.
Spherical and ellipsoidal (including flat ellipsoidal) uncertainty sets were considered in









[13]-[15], [44], [45]. The robust Capon beamforming approaches presented in [44], [45]
coupled the spatial filtering formulation of the standard Capon beamformer (SCB)

in [10] with a spherical or ellipsoidal uncertainty set of the array steering vector.

We coupled the covariance fitting formulation of SCB in [12] with an ellipsoidal or

spherical uncertainty set, to obtain a robust Capon beamformer (RCB) [13], [14].

Interestingly, the methods in [13], [14], [44], [45] turn out to be equivalent and

to belong to the extended class of diagonal loading approaches, but the correspond-

ing amount of diagonal loading can be calculated precisely based on the ellipsoidal

uncertainty set of the array steering vector. However, our RCB in [13] is simpler

and computationally more efficient than its equivalent counterparts; and its compu-

tational complexity is comparable to that of SCB. Moreover, our RCB gives a simple

way of eliminating scaling ambiguity, when estimating the power of the desired signal.

The approaches in [44], [45] did not consider the scaling ambiguity problem.

1.3 Scope of the Work

Our study focused on theory and applications of the RCB algorithm, which

is a natural extension of the covariance fitting formulation of SCB to the case of

uncertain array steering vectors. We derived RCB for two cases: non-degenerate

ellipsoidal constraints on the steering vector; and flat ellipsoidal constraints. We pre-

sented a diagonal loading interpretation of RCB, and showed that its diagonal loading

level can be calculated precisely based on the uncertainty set of the steering vector.

We also provided insights into relationships among the recent three robust adaptive

beamformers. In addition, we address several extensions of RCB for aero-acoustic

and ultrasonic applications. We offer numerous simulated and experimental examples

to demonstrate the efficacy of RCB and also its extensions. More applications and

analyses on RCB can be found in [63]-[66].














CHAPTER 2
ROBUST CAPON BEAMFORMING (RCB)

In this chapter, we first formulate the problem of interest. Then we briefly

review two formulations of the standard Capon beamformer, namely the spatial fil-

tering SCB and the covariance fitting SCB, and demonstrate their equivalence. Next,

we focus on the robust Capon beamformer by considering two cases( i.e., the case of

non-degenerate ellipsoidal constraints on the steering vector and the case of flat ellip-

soidal constraints). These two cases are treated separately because of the differences

in their detailed computational steps and in the possible values of the associated

Lagrange multipliers. We will also provide a diagonal loading interpretation of RCB.

2.1 Problem Formulation

Consider an array comprising M sensors. Let R denote the theoretical covari-

ance matrix of the array output vector. We assume that R > 0 (positive definite)

has the following form:
K
R = a0aoa0 + aaka* + Q (2.1)
k=1

where (a2, {a}kK=i) are the powers of the (K + 1) uncorrelated signals impinging

on the array, (ao, {ak}K=1) are the so-called steering vectors that are functions of the

location parameters of the sources emitting the signals (e.g., their directions of arrival

(DOAs)), (.)* denotes the conjugate transpose, and Q is the noise covariance matrix

(the "noise" comprises nondirectional signals; and hence Q usually has full rank as

opposed to the other terms in (2.1) whose rank is equal to one). We assume that the

first term in (2.1) corresponds to the signal-of-interest (SOI), and that the remaining









rank-one terms correspond to K interference. We assume that


lao 12 = M (2.2)

where [| | denotes the Euclidean norm. We note that the above expression for R

holds for both narrowband and wideband signals. For narrowband signals, R is the

covariance matrix at the center frequency. For wideband signals, R is the covariance

matrix at the center of a given frequency bin. Let

R = uru* (2.3)


where the columns of U contain the eigenvectors of R; and the diagonal elements of

the diagonal matrix r, 71 >Y 72 > ... 7M, are the corresponding eigenvalues. In

practical applications, R is replaced by the sample covariance matrix R, where
N
R= EYnY (2.4)
n=l
N is the number of snapshots and Yn is an M x 1 vector representing the nth snapshot

with the form

Yn = aoso(n) + en, (2.5)

with so(n) denoting the waveform of the SOI, and en being the interference-plus-

noise vector for the nth snapshot. In the radar application, the SOI-free data vectors

are available, so the interference-plus-noise covariance matrix Rin may be estimated

and used for beamforming. However, in many other applications, the SOI is always

present in the data. Hence we cannot exploit the interference-plus-noise covariance

matrix. Instead, we design beamformers based on R. Throughout our study, we

considered this more general circumstance.

The robust adaptive beamforming problem we will deal with in this chapter

is as follows: extend the Capon beamformer so as to be able to accurately determine

the power of SOI, even when knowledge of its steering vector, ao, is imprecise. More









specifically, we assume that our only knowledge about ao is that it belongs to an

uncertainty ellipsoid. The non-degenerate ellipsoidal uncertainty set has the following

form

a I (a- a)*C- a- ) 1 (2.6)

where a (the assumed steering vector of SOI) and C (a positive definite matrix) are

given. The ellipsoidal uncertainty set above is equivalent to

{a a=C u+ a, Iul| < 1 (2.7)

where CA is the Hermitian square root matrix of C (i.e., C2(C )* = C). If we
assume that u is a circularly symmetric random vector, with a mean of zero and

covariance matrix equal to identity matrix I, then we can prove the following results:

E [a] = a (2.8)

and

E [(a a)(a a)*] = C (2.9)

The ellipsoidal uncertainty set may be generated from the mean and covariance ma-
trix of the array steering vector, calculated based on the measurements from repeated
trials in the experiments. Even without the above statistical assumptions, we can

still exploit semidefinite programming to obtain a tight ellipsoidal uncertainty set,
the minimum volume ellipsoid covering all the samples of the array steering vectors

[45]. The case of a flat ellipsoidal uncertainty set is considered in Section 2.3.2.
In particular, if C in (2.6) is a scaled identity matrix, (i.e., C = el), we have
the following uncertainty sphere:

{a I Ila all2 } (2.10)

where E is a user parameter that can be calculated based on a priori knowledge of the
array steering vector. For example, if we know the varying ranges of the gain, phase,









and position errors of each sensor, we can determine the corresponding uncertainty

region, which is an annulus sector. To illustrate this point, we use a uniform linear

array. Assume that at the mth sensor, the uncertainty is caused by gain error rm,

phase error 7m, and element position error Am and Aym. Thus we have

am = dm(1 + rm)eipm, m = 1,..., M (2.11)

where the total phase error Pm is

Axm sin 0 + Aym cos 0
Pm = 7m + Wc (2.12)
c

with 0 being the signal direction of arrival, wc representing the carrier frequency,

and c denoting the wave propagating velocity. Assume that the varying ranges for

rm, 7m, Axm, Aym are known. Then we can determine the uncertainty region Qm for

am and the maximum total error em, which is given by

em = max am amI, m= 1, M (2.13)
amEnm

Hence we can choose e as follows:
M
E- Me12 (2.14)
m=1
From an application standpoint, RCB makes it possible to avoid the time-consuming

array calibration process, provided that the steering vector uncertainty set is known.

Our study focused on estimating the SOI power ao from R (or more practically

R) when the knowledge of ao is imprecise. However, our beamforming approaches

we present herein can also be used for other applications, such as signal waveform

estimation [44], [45].

2.2 Standard Capon Beamforming
2.2.1 Spatial Filtering SCB

The common formulation of the beamforming problem that leads to the spatial

filtering form of SCB [5], [10], [11] is as follows.









(a) Determine the M x 1 weight vector wo that solves the following linearly

constrained quadratic problem:

minw*Rw subject to w*ao = 1 (2.15)
w

(b) Use w*Rwo as an estimate of a2.

The solution to (2.15) is easily derived:

R-lao
wo = R-l (2.16)
a*R-lao
Using (2.16) in Step (b) yields the following estimate:

a2 = a0 (2.17)
0 aOR-lao
Note that (2.15) can be interpreted as an adaptive spatial filtering problem:

given R and ao, we wish to determine the weight vector wo as a spatial filter that can

pass SOI without distortion; and at the same time minimize undesirable interference

and noise contributions in R.

2.2.2 Covariance Fitting SCB

The Capon beamforming problem can also be reformulated into a covariance

fitting form. To describe the details of our approach, we first proved that a2 in (2.17)

is the solution to the following problem [12], [14]:

maxa2 subject to R- a2aoa O > 0 (2.18)

where the notation A > 0 (for any Hermitian matrix A) means that A is positive

semi-definite. The previous claim follows from the following readily verified equiva-

lences (here R-1/2 is the Hermitian square root of R-'):

R a02aa* > 0

I a2R-1/2aoaR-12 > 0 ::

1 a2aR-lao > 0 <>

o2 < aR a = 2 (2.19)
Sa;R-la









Hence a2 = -2 is indeed the largest value of a2 that satisfies the constraint in (2.18).

Note that (2.18) can be interpreted as a covariance fitting problem: given R

and ao we wish to determine the largest possible SOI term (a2aoa*) that can be a part

of R, under the natural constraint that the residual covariance matrix be positive

semi-definite.

2.3 Robust Capon Beamforming

The robust Capon beamformer is derived by a natural extension of the covari-

ance fitting SCB in Section 2.2.2 to the case of uncertain steering vector. In doing so

we directly obtain a robust estimate of a2, without any intermediate calculation of a

vector w [13], [14].

2.3.1 Non-Degenerate Ellipsoidal Uncertainty Set

When the uncertainty set of the steering vector a is a non-degenerate ellipsoid

as in (2.6), the RCB problem has the following form [13], [14]:

max a2 subject to R a2aa* > 0
Cr2,a
for any a satisfying (a a)* C-1 (a a) < 1 (2.20)


where a and C are given.

The RCB problem in (2.20) can be readily reformulated as a semi-definite

program (SDP) [14]. Indeed, using a new variable X = 1/a2 along with the standard

technique of Schur complements [5], [67] we can rewrite (2.20) as:

min X subject to
xa

R a
a* x -

[( aC ) a-]a]o (2.21)
( a) TY 1 >0 (2.21)

The constraints in (2.21) are so-called linear matrix inequalities, and hence

(2.21) is an SDP, which requires O(oM6) flops if the SeDuMi type of software [68]









is used to solve it, where g is the number of iterations, usually on the order of v/iM.

However, the approach we present below only requires 0(M3) flops.

For any given a, the solution &- to (2.20) is indeed given by the counterpart of

(2.17) with ao replaced by a, as shown in Section 2.2.2. Hence (2.20) can be reduced

to the following problem

mina*R-la subject to (a a)* C-1 (a a) < 1 (2.22)
a

To exclude the trivial solution a = 0 to (2.20), we assume that

a*C-a > 1 (2.23)

Note that we can decompose any matrix C > 0 in the form:

1
C-1= 1D*D (2.24)
E

where for some e > 0,

D = vC-1/2 (2.25)

Let

S= Da, =Da, R =DRD* (2.26)

Then (2.22) becomes

mina*lR-i subject to a 2 < e (2.27)

Hence without loss of generality, we will consider solving (2.22) for C = el, i.e.,

solving the following quadratic optimization problem under a spherical constraint:

mina*R-la subject to Ila a112 < e (2.28)
a

To exclude the trivial solution a = 0 to (2.28), we now need to assume that


l|all2 > e


(2.29)








Let S be the set defined by the constraints in (2.28). To determine the solution
to (2.28) under (2.29), consider the function:

hi(a, v) = a*R-la + v (la a|12 e) (2.30)

where v > 0 is the real-valued Lagrange multiplier satisfying R-' + vI > 0 so
that the above function can be minimized with respect to a. Evidently we have
hi(a, v) < a*R-la for any a E S with equality on the boundary of S. Equation
(2.30) can be written as

hi (a, v)= a 1) (R- + VI) a +(R a

-v2 a*R-1 + vI)-a + va*a ve (2.31)

Hence the unconstrained minimization of h (a, v) w.r.t. a, for fixed v, is given by
(R-1 -1
o = + a (2.32)

= (I + R)- a (2.33)

where we have used the matrix inversion lemma [5] to obtain the second equality.
Clearly, we have

h2(v) h(, v) (R1 I v (2.34)

< a*R-a for any a E S (2.35)

Maximization of h2(v) with respect to v gives

(V- +l (2.36)

which indeed satisfies

Il|o a1|2 = (2.37)

Hence do belongs to the boundary of S. Therefore, ao is the sought solution.








Using (2.33) in (2.37), the Lagrange multiplier v > 0 is then obtained as the
solution to the constraint equation:

h2(v) I(I R)- -12 = (2.38)

Making use of R = UPU* and z = U*a, (2.38) can be written as

h2m (1 )2 (2.39)
m=l(1 + Vm)
Note that h2(v) is a monotonically decreasing function of v > 0. According to (2.29)
and (2.38), h2(0) > E and hence v 0 0. From (2.39), it is clear that lim,,,o h2(v) =
0 < e. Hence there is a unique solution v > 0 to (2.39). By replacing the ym in (2.39)
with YM and -y, respectively, we can obtain the following tighter upper and lower
bounds on the solution v > 0 to (2.39):

all- < v < (2.40)
71C 7YM V
By dropping the 1 in the denominator of (2.39), we can obtain another upper bound
on the solution v to (2.39):
1
v< 1im (2.41)

The upper bound in (2.41) is usually tighter than the upper bound in (2.40) but not
always. In summary, the solution v > 0 to (2.39) is unique and it belongs to the
following interval:

Hal-<- < v < mmin -M IZ 2 Ha V- (2.42)
( 1=1
Once the Lagrange multiplier v is determined, ao is determined by using (2.33)
and &02 is computed by using (2.17) with ao replaced by Ao. Hence the major compu-
tational demand of our RCB comes from the eigen-decomposition of the Hermitian
matrix R, which requires O(M3) flops. Therefore, the computational complexity of
our RCB is comparable to that of the SCB.









Next observe that both the power and the steering vector of SOI are treated
as unknowns in our robust Capon beamforming formulation (see (2.20)), and hence

that there is a "scaling ambiguity" in the SOI covariance term in the sense that

(a2, a) and (a2/a, a1/2a) (for any a > 0) give the same term a2aa*. To eliminate
this ambiguity, we use the knowledge that l|ao 12 = M (see (2.2)) and hence estimate

o'2 as [14]
-2
U2o = &jllaoll2/M (2.43)

where &2 is obtained via replacing ao in (2.17) by do in (2.32). The numerical exam-
2
ples in [14] confirm that 0o is a (much) more accurate estimate of ao than -2. To

summarize, our proposed RCB approach consists of the following steps.

The Proposed RCB (spherical constraint)

Step 1: Compute the eigen-decomposition of R (or more practically of R).

Step 2: Solve (2.39) for v, e.g., by a Newton's method, using the knowledge that

the solution is unique and it belongs to the interval in (2.42).

Step 3: Use the v obtained in Step 2 to get

ao = a U (I + vr)-' U*a (2.44)

where the inverse of the diagonal matrix I + vF is easily computed. (Note that

(2.44) is obtained from (2.33).)

Step 4: Compute &2 by using

&0 = 1 (2.45)
o *Ur (v-2I + 2v-Fl + r2)-1 U*a

where the inverse of v-2I + 2v-lr + r2 is also easily computed. Note that ao

in (2.17) is replaced by do in (2.32) to obtain (2.45). Then use the &2 in (2.43)
S
to obtain o0 as the estimate of ao.









We remark that in all of the steps above, we do not need to have 7m > 0 for
all m = 1, 2, ... M. Hence R or R can be singular, which means that we can allow
N < M to compute R.
In other applications, such as communications, the focus is on SOI waveform
estimation. Let so(n) denote the waveform of the SOI. Then once we have estimated
the SOI steering vector with our RCB, so(n) can be estimated like in the SCB as
follows:

o(n) = w x, (2.46)

where ao in (2.32) is used to replace ao in (2.16) to obtain Wro:

R-lfo
w0 =- R (2.47)

(R + 11) A
(= R+ )1a (2.48)
d* (R + 4 )-1 R (R + -I) a

Note that our robust Capon weight vector has the form of diagonal loading except for
the real-valued scaling factor in the denominator of (2.48). However, the scaling factor
is not really important since the quality of the SOI waveform estimate is typically
expressed by the signal-to-interference-plus-noise ratio (SINR)

a w21**ao12
SINR = ao a12 (2.49)
w (k=1 Ukakaa + Q) wo

which is independent of the scaling of the weight vector.
When C is not a scaled identity matrix, the diagonal loading is added to the
weighted matrix R defined in (2.26) and we refer to this case as extended diagonal
loading. To exclude the trivial solution a = 0 to (2.20), we now need to assume, like
in (2.29), that

11112 > E (2.50)

which is equivalent to
*C-li^ > 1 (2.51)









The discussions above indicate that our robust Capon beamforming approach

belongs to the class of (extended) diagonally loaded Capon beamforming approaches.

However, unlike earlier approaches, our approach can be used to determine exactly

the optimal amount of diagonal loading needed for a given ellipsoidal uncertainty set

of the steering vector, at a very modest computational cost.

Our approach is different from the recent RCB approaches in [44], [45]. The

latter approaches extended Step (a) of the spatial filtering SCB in Section 2.2.1 to

take into account the fact that when there is uncertainty in ao, the constraint on

w*ao in (2.15) should be replaced with a constraint on w*a for any vector a in the

uncertainty set (the constraints on w*a used in [44] and [45] are different from one

another); then the so-obtained w is used in w*Rw to derive an estimate of ao, as

in Step (b) of the spatial filtering SCB. Unlike our approach, the approaches of [44]

and [45] do not provide any direct estimate 0o. Hence they do not provide a simple

way (such as (2.43)) to eliminate the scaling ambiguity of the SOI power estimation

that is likely a problem for all robust beamforming approaches (this problem was in

fact ignored in both [44] and [45]). Yet SOI power estimation is often the main goal

in many applications including radar, sonar, acoustic and medical imaging.
Despite the apparent differences in formulation, we prove in Appendices A

and C that our RCB gives the same weight vector as the approaches presented in

[44], [45], yet our RCB is computationally more efficient [13]. The approach in [44]
requires O(QM3) flops [69], where Q is the number of iterations, usually on the order

of V/M, whereas our RCB approach requires O(M3) flops. Moreover, our RCB can

be readily modified for recursive implementation by adding a new snapshot to R

and possibly deleting an old one. By using a recursive eigen-decomposition updating

method [70], [71] with our RCB, we can update the power and waveform estimates in

O(M2) flops. No results are available so far for efficiently updating the second-order
cone program (SOCP) approach in [44]. The approach in [45] can be implemented









recursively by updating the eigen-decomposition similarly to our RCB. However, its

total computational burden can be higher than for ours, as explained in the next

subsection.

We also show in Appendix B that, although this aspect was ignored in [44],

[45], the approaches presented in [44], [45] can also be modified to eliminate the
scaling ambiguity problem that occurs when estimating the SOI power a2 [13].

2.3.2 Flat Ellipsoidal Uncertainty Set

When the uncertainty set of a is a flat ellipsoid, as is considered in [45] to

make the uncertainty set as tight as possible (assuming that the available a priori

information allows that), (2.20) becomes [13]

max a2 subject to R a2aa* > 0
02,a
a=Bu + lu < 1 (2.52)

where B is an M x L matrix (L < M) with full column rank and u is an L x 1 vector.

(When L = M, the second constraint in (2.52) becomes (2.6) with C = BB*.) Below

we provide a separate treatment of the case of L < M due to the differences from

the case of L = M in the possible values of the Lagrange multipliers and the detailed

computational steps. The RCB optimization problem in (2.52) can be reduced to

(see (2.22)):
min(Bu + a)*R-1(Bu + a) subject to Ilull < 1 (2.53)
U
Note that

(Bu + a)*R-l(Bu + a) = u*B*R-1Bu + a*R-1Bu + u*B*R-1a + a*R-'a (2.54)

Let

R = B*R-1B > 0 (2.55)

and


(2.56)









Using (2.54)-(2.56) in (2.53) gives

minu*Ru + *u + u* subject to |luI < 1 (2.57)
U

To avoid the trivial solution a = 0 to the RCB problem in (2.52), we im-
pose the following condition (assuming fi below exists, otherwise there is no trivial
solution). Let fi be the solution to the equation

Bfi a = 0 (2.58)

Hence
i = -Bta (2.59)

Then we require that
i*Bt*Bti > 1 (2.60)

where Bt denotes the Moore-Penrose pseudo-inverse of B.
The Lagrange multiplier methodology can again be used to solve (2.52) [72].
Let

hi (u, 0) = u*Ru + A*u + u*a + P(u*u 1) (2.61)

where i > 0 is the Lagrange multiplier [73]. Differentiation of (2.61) with respect to
u gives
R~i + A + oifi = 0 (2.62)

which yields
fi = -(R + I)- A (2.63)

If IlR-1all < 1, then the unique solution in (2.63) with 0 = 0, which is it = -R-l',
solves (2.57). If Ili-1'AI > 1, then f > 0 is determined by solving

2( I (ft P-I)1 2 = 1 (2.64)

Note that h2(0) is a monotonically decreasing function of 0 > 0. Let


(2.65)


R = Ot"*








where the columns of U contain the eigenvectors of R and the diagonal elements of
the diagonal matrix f, j1 > 2 > L, are the corresponding eigenvalues. Let

z = T0* (2.66)

and let it denote the Ith element of 2. Then

L 112
2 + Zp 2 1 (2.67)
1=1

Note that lim,,oo2(i) = 0 and h2(0) = I|Ri-l1I > 1. Hence there is a unique
solution to (2.67) between 0 and oo. By replacing the 7t in (2.67) with jL and 71,
respectively, we obtain tighter upper and lower bounds on the solution to (2.67):

I1a| a1 111 (2.68)

Hence the solution to (2.67) can be efficiently determined by using, for example, the
Newton's method, in the above interval. Then the solution 0 to (2.67) is used in
(2.63) to obtain the fi that solves (2.57).
To summarize, our proposed RCB approach consists of the following steps.
The Proposed RCB (flat ellipsoidal constraint)

Step 1: Compute the inverse of R (or more practically of R) and calculate R and
a using (2.55) and (2.56), respectively.

Step 2: Compute the eigen-decomposition of R (see (2.65)).

Step 3: If IIR-1ll < 1, then set 0 = 0. If 11R-111 > 1, then solve (2.67) for 0, e.g.,
by a Newton's method, using the knowledge that the solution is unique and it
belongs to the interval in (2.68).

Step 4: Use the v obtained in Step 3 to get


f = -U (f + I)-10* A


(2.69)









(which is obtained from (2.63)). Then use the fi to obtain the optimal solution

to (2.52) as

ao = Bf + a (2.70)

Step 5: Compute 802 by using (2.17) with ao replaced by ao and then use the &2 in

(2.43) to obtain the estimate of o02.

Hence, under the flat ellipsoidal constraint the complexity of our RCB is also

O(M3) flops, which is on the same order as for SCB and is mainly due to computing

R-1 and the eigen-decomposition of R. If L << M, then the complexity is mainly

due to computing R-1. Note, however, that to compute i, we need O(L3) flops while

the approach in [45] requires 0(M3) flops (and L < M).

2.4 Diagonal Loading Interpretation of RCB

In many applications, such as in communications or the global positioning

system, the focus is on SOI waveform estimation. The waveform of the SOI, so(n),

as in (2.5) can be estimated as follows:

o(n) = *Yn (2.71)

where is the corresponding weight vector. For RCB, we can substitute the esti-

mated steering vector ao in lieu of ao in (2.16) to obtain w.

Diagonal loading is a popular approach to mitigate the performance degrada-

tions of SCB in the presence of steering vector error or the small sample size problem.

As the name implies, its weight vector has a diagonally loaded form:

w = ,(R + SI)- a (2.72)

where 6 denotes the diagonal loading level. Also, in (2.72) K is a scaling factor,

which can be important for accurate power estimation; however, it is immaterial

for waveform estimation since the quality of the SOI waveform estimate is typically









measured by the signal-to-interference-plus-noise ratio (SINR)

SINR = ,0l2-*ao12 (2.73)
(E= oKaka + Q)
which is independent of i.
As a matter of fact, RCB can be interpreted in the unified framework of
diagonal loading based approaches.

2.4.1 Nondegenerate Ellipsoidal Uncertainty Set

Using ao in (2.32) to replace ao in (2.16), we can obtain the following RCB
weight vector for the case of non-degenerate ellipsoidal uncertainty set:

S(R + I)-1 -
Wad = (R + 4I)- R (R + -I)1 (

When C is not a scaled identity matrix, the diagonal loading is added to the
weighted matrix 1i defined in (2.26) instead of R and we refer to this case as the
extended diagonal loading.

2.4.2 Flat Ellipsoidal Uncertainty Set

The RCB weight vector for the case of flat ellipsoidal Uncertainty set has the
form:

-1 0
WRCB_F --1-
0adR-ldo
(R + 1BB*)- a
= (2.75)
a* (R + BB*)-R (R + BB*)-
To obtain (2.75) we have used the fact (also using (2.63) in (2.70)) that

R- o = -R-1B(R+ PI)- + R-la

= -R-B(B*R-1B + PI)- B*R-1a + R-a

= (R+ BB* a (2.76)

where the last equality follows from the matrix inversion lemma. We see that in this
case, the RCB weight vector again has an extended diagonally loaded form.









Despite the differences in the formulation of our RCB problem and that in [45],

we can prove that the WRCB.F in (2.75) and the optimal weight in [45] are identical

(Appendix C).

2.5 Numerical Examples

Next, we provide numerical examples to compare the performances of the

SCB and RCB. In all of the examples considered below, we assume a uniform linear

array with M = 10 sensors and half-wavelength sensor spacing, and a spatially white

Gaussian noise whose covariance matrix is given by Q = I.

Example 2.1: Comparison of SCB and RCB for the case of finite

number of snapshots without look direction errors. We consider the effect of

the number of snapshots N on the SOI power estimate when the sample covariance

matrix R in (2.4) is used in lieu of the theoretical array covariance matrix R in both

the SCB and RCB. (Whenever R is used instead of R, the average power estimates

from 100 Monte-Carlo simulations are given. However, the beampatterns shown

are obtained using R from one Monte-Carlo realization only.) The power of SOI is
4a = 10 dB and the powers of the two (K = 2) interference assumed to be present

are a2 = a' = 20 dB. We assume that the steering vector uncertainty is due to the

uncertainty in the SOI's direction of arrival o0, which we assume to be 0o + A. We

assume that a(o0) belongs to the uncertainty set


I|a(0o) a|l2 < e; = a(0o + A) (2.77)

where e is a user parameter. Let eo = Ila(Oo) a112. Then choosing e = co gives the

smallest set that includes a(0o). However, since A is unknown in practice, the e we
choose may be greater or less than co. To show that the choice of e is not a critical

issue for our RCB approach, we will present numerical results for several values of e.

We assume that the SOI's direction of arrival is 0o = 0 and the directions of arrival

of the interference are 01 = 600 and 02 = 80.







28



18 o 8- o 0
A
S6 S6-
4- 4-
E 2- 2-
t 0- o 0 o0
0 -2- 0 -2-

S -4- ( -4-
-6 A RCB (Sample R) -6 A RCB (Sample R)
o SCB (Sample R) o SCB (Sample R)
1- SCB (Theoretical R) -SCB (Theoretical R)
-10 I I I I I I -110 I I I I I I
10 20 30 40 50 60 70 80 90 100 10 20 30 40 50 60 70 80 90 100
Number of Snapshots Number of Snapshots
(a) (b)

2 2

versus N. (a) e = 0.5. (b) e = 3.5. The true SOI power is 10 dB and Eo = 0 (i.e., no
mismatch).

2 2
In Figure 2.1, we show uo2 and o-o versus the number of snapshots N for the

no mismatch case; hence A = 0 in (2.77) and consequently eo = 0. Note that the

power estimates obtained by using R approach those computed via R as N increases,

and that our RCB converges much faster than the SCB. The SCB requires that N

is greater than or equal to the number of array sensors M = 10. However, our RCB

works well even when N is as small as N = 2.

Figure 2.2 shows the beampatterns of the SCB and RCB using R as well as

R with N = 10, 100, and 8000 for the same cae as in Figure 2.1. Note that the

weight vectors used to calculate the beampatterns of RCB in this example (as well

as in the following examples) are obtained by using the scaled estimate of the array

steering vector /Mco/nsoj|l in (2.16) instead of Ao. The vertical dotted lines in the

figure denote the directions of arrival of the SOI and the interference. The horizontal

dotted lines in the figure correspond to 0 dB. Figure 2.2(a) shows that although the

RCB beampatterns do not have nulls at the directions of arrival of the interference as

deep as those of the SCB, the interference (whose powers are 20 dB) are sufficiently

suppressed by the RCB to not disturb the SOI power estimation. Regarding the poor








29
















Si 40 II II I 4 0 I

20- 20- "..-.. .. -

I -. ;. '"
0 -- ...... ..... .. ....... .......... ........B ; o .... ........ ........... ....... ,--.R. ... .
E E
... .. i i . . i I i











(a) (b)
-20 -20

-40 -40 -
-
-60- -60
S-80 80

-100 S1B I I I I
-80 -60 -40 -20 0 20 40 60 80 -80 -60 -40 -20 0 20 40 60 80
0 degree degree

(a) (b)
S 0 dB ad 40 (o

20- 20



EE


S-40 -40

-8o -0 -0

-80 -80-
0I I I I -100 I I I I
-80 -60 -40 -20 0 20 40 60 80 -80 -60 -40 -20 0 20 40 60 80
degree degree

(c) (d)

Figure 2.2: Comparison of the beampatterns of SCB and RCB when E = 3.5. (a)
Using R. (b) Using R with N = 10. (c) Using R with N = 100. (d) Using R with
N = 8000. The true SOI power is 10 dB and eo = 0 (no mismatch).











,.. I I .

8 8-
8 -A8B
6 S 6-

E
2-.
_0 - o o o o o-o-o o- -o- .... - -- - -
I .2-
4 -4- -4-
-6 A RCB (Sample R) -6 A RCB (Sample R)
o SCB (Sample R) o SCB (Sample R)
-8 -- RCB (Theoretical R) -8 RCB (Theoretical R)
S--- SCB (Theoretical R) SCB (Theoretical R)
-1 C -1 C
10 20 30 40 50 60 70 80 90 100 10 20 30 40 50 60 70 80 90 100
Number of Snapshots Number of Snapshots
(a) (b)

Figure 2.3: Comparison of 02 (SCB using R and R) and -o (RCB using R and R)
versus N. (a) e = 2.5. (b) e = 4.5. The true SOI power is 10 dB and Co = 3.2460
(corresponding to A = 2.00).


performance of SCB for small N, note that the error between R and R can be viewed

as due to a steering vector error [38].

Example 2.2: Comparison of SCB and RCB for the case of finite

number of snapshots in the presence of look direction errors. This example

is similar to Example 2.1 except that now the mismatch is A = 20 and accordingly

Eo = 3.2460. We note from Figure 2.3 that even a relatively small A can cause a

significant degradation of the SCB performance. As can be seen from Figure 2.4, the

SOI is considered to be an interference by SCB and hence it is suppressed. On the
,2
other hand, the SOI is preserved by our RCB and the performance of o-o obtained

via our approach is quite good for a wide range of values of e. Note that the RCB

also has a smaller "noise gain" than the SCB.

Example 2.3: Comparison of the RCB method and a fixed diagonal

loading level based approach. In Figure 2.5, we compare the performance of

our RCB with a fixed diagonal loading level based approach. Specifically, the fixed

loading level was chosen equal to 10 times the noise power (assuming the knowledge

of the noise power). Consider the same case as Figure 2.4(d) except that now we













































) -60 -30 0
0 degree


30 60 9


C





ca

|-6C


-8C


-1-0 I I --


S -60 -30 0
0 degree

(b)


0 30 60 90 -'90 -60 -30 0
0 degree 0 degree

(c) (d)


Figure 2.4: Comparison of the beampatterns of SCB and RCB when e = 1.0 when

(a) using R and (b) using R with N = 10, and when e = 4.5 for (c) using R and (d)

using R with N = 10. The true SOI power is 10 dB and co = 3.2460 (corresponding

to A = 2.00).


_ -- RCB
S i



I I I


. ......... .. .... .. .. ............ ...


i









|-- RCB |
I -- SCB .


30 60


2n ..i i


) -9









assume that R is available and we vary the SNR by changing the SOI or noise power.

For Figures 2.5(a), 2.5(c) and 2.5(e), we fix the noise power at 0 dB and vary the SOI

power between -10 dB and 20 dB. For Figures 2.5(b), 2.5(d) and 2.5(f), we fix the

SOI power at 10 dB and vary the noise power between -10 dB and 20 dB. Figures

2.5(a) and 2.5(b) show the diagonal loading levels of our RCB as functions of the

SNR. Figures 2.5(c) and 2.5(d) show the SINRs of our RCB and the fixed diagonal

loading level approach and Figures 2.5(e) and 2.5(f) show the corresponding SOI

power estimates, all as functions of the SNR. Note from Figures 2.5(a) and 2.5(b)

that our RCB adjusts the diagonal loading level adaptively as the SNR changes. It

is obvious from Figure 2.5 that our RCB significantly outperforms the fixed diagonal

loading level approach when the SNR is medium or high.

Example 2.4: Comparison of RCB, SCB and the delay-and-sum

method in the presence of array calibration errors. We consider an imaging

example, where we wish to determine the incident signal power as a function of the

steering direction 0. We assume that there are five incident signals with powers 30,

15, 40, 35, and 20 dB from directions -35, -15, 0, 100, and 400, respectively.

To simulate the array calibration error, each element of the steering vector for each

incident signal is perturbed with a zero-mean circularly symmetric complex Gaussian

random variable so that the squared Euclidean norm of the difference between the

true steering vector and the assumed one is 0.05. The perturbing Gaussian random

variables are independent of each other.

Figure 2.6 shows the power estimates of SCB and RCB, obtained using R, as

a function of the direction angle, for several values of e. The small circles denote the

true (direction of arrival, power)-coordinates of the five incident signals. Figure 2.6

also shows the power estimates obtained with the data-independent beamformer using

the assumed array steering vector as the weight vector. This approach is referred to

as the delay-and-sum beamformer. We note that SCB can still give good direction of








































SNR (dB) SNR

(a) (b)


-10 -5 0 5
SNR (dB)


5
SNR (dB)

(e)


10 15 20


A



0 0
o o
S10 -

Z
5-
0
a

0-



-0 -5 0 5 101 5 20
SNR (dB)

(d)
15 I I I



0



0-



0
6 0



CI)
O o

-5-
o Fixed diagonal loading
A RCB
-0 SOI power


20 -10 -5 0 5
SNR (dB)

(f)


10 15 20


Figure 2.5: Comparison of a fixed diagonal loading level approach and our RCB when

e = 4.5 and eo = 3.2460 (corresponding to A = 2.00). (a) (c) (e) Signal power change.

(b) (d) (f) Noise power change.


o Fixed diagonal loading
A RCB


o Fixed diagonal loading


A

A
a RCB&





A
0 0
0

'0
S0
a
6
o


I I










RCB -- RCB
40- SCB 40-, -- SCB
S- Delay-and-sum - Delay-and-sum


E 2 i s
S





-10 -10
CL a



-60 -40 -20 0 20 40 60 -60 -40 -20 0 20 40 60
0 degree 0 degree
(a) (b)

Figure 2.6: Power estimates (using R) versus the steering direction 0. (a) e = 0.03.
(b) e = 0.1. The true powers of the incident signals from -350, -150, 0, 100, and
400 are denoted by circles, and co = 0.05.


arrival estimates for the incident signals based on the peak power locations. However,

the SCB estimates of the incident signal powers are way off. On the other hand, our

RCB provides excellent power estimates of the incident sources and can also be used

to determine their directions of arrival based on the peak locations. The delay-and-

sum beamformer, however, has much poorer resolution than both SCB and RCB.

Moreover, the sidelobes of the former give false peaks.

Example 2.5: Comparison of SCB, RCB with spherical constraint

and RCB with flat ellipsoidal constraint in the presence of look direction

errors. We examine now the effects of the spherical and flat ellipsoidal constraints on

SOI power estimation. We consider SOI power estimation in the presence of several

strong interference. We will vary the number of interference from K = 1 to K = 8.

The power of SOI is a2 = 20 dB and the interference powers are of = * = a2 = 40

dB. The SOI and interference directions of arrival are 0o = 10, 01 = -75, 82 =

-60, 03 = -450, 4 = -300, 05 = -100, 06 = 25, 07 = 35, 08 = 500. We assume

that there is a look direction mismatch corresponding to A = 20 and accordingly

eo = 3.1349.













2C-:::l:::^A--- -----A----- ------- ------ -~ 2 ----s=.---ii ---r----a- ----&---~ -

S 20


10 310
0



ri---"T ---------'--- l--- l --t--- " "' I -- --
15 1


O 6



1 2 3 4 5 6 7 8 1 2 3 4 5 6 8
Number of Interferences Number of Interferences
(a) (b)

Figure 2.7: Comparison of 0a (SCB), Oo (RCB with flat ellipsoidal constraint with
2
L = 2), and a0 (RCB with spherical constraint), based on R, versus the number
of interference K. (a) 6 = 1.8. (b) 6 = 2.40. The true SOI power is 20 dB and
co = 3.1349 (corresponding to A = 2).


Figure 2.7 shows the SOI power estimates, as a function of the number of

interference K, obtained by using SCB, RCB (with flat ellipsoidal constraint), and

the more conservative RCB (with spherical constraint) all based on the theoretical

array covariance matrix R. For RCB with flat ellipsoidal constraint, we let B contain

two columns with the first column being a(0o + A) a(0o + A 6) and the second

column being a(0o + A) a(00 + A + 3). Note that choosing 6 = A = 2' gives the

smallest flat ellipsoid that this B can offer to include a(00). However, we do not know

the exact look direction mismatch in practice. We choose 6 = 1.80 and 6 = 2.40 in

Figures 2.7(a) and (b), respectively. For RCB with spherical constraint, we choose e

to be the larger of Ila(o +A) -a(0 +A-3)112 and Ila(o +A) -a(0 +A+ )11 2. Note

that RCB with flat ellipsoidal constraint and RCB with spherical constraint perform

similarly when K is small. However, the former is more accurate than the latter for

large K.

Figure 2.8 gives the beampatterns of the SCB and RCBs using R as well

as R with N = 10 for various K. For large K, the more conservative RCB with

spherical constraint amplifies the SOI while attempting to suppress the interference,


-25 A RB (Flat ellipsoid)
-- RCB (Sphere)
-o- SCB _-a-- ---
..... ...----------


-A R B (Flat ellipsoid)
-. RCB (Sphere)
-- SCB --.----
1 -------- -- ---- ---- ---'









as shown in Figure 2.8. On the other hand, the RCB with flat ellipsoidal constraint

maintains an approximate unity gain for the SOI and provides much deeper nulls for

the interference than the RCB with spherical constraint at a cost of worse noise gain.

As compared to the RCBs, the SCB performs poorly as it attempts to suppress the

SOI. Comparing Figures 2.8(b) with 2.8(a), we note that for small K and N, RCB

with spherical constraint has a much better noise gain than RCB with flat ellipsoidal

constraint, which has a better noise gain than SCB. From Figure 2.8(d), we note that

for large K and small N, RCB with flat ellipsoidal constraint places deeper nulls at

the interference angles than the more conservative RCB with spherical constraint.

Figure 2.9 shows the SOI power estimates versus the number of snapshots

N for K = 1 and K = 8 when the sample covariance matrix R is used in the

beamformers. Note that for small K, RCB with spherical constraint converges faster

than RCB with flat ellipsoidal constraint as N increases, while the latter converges

faster than SCB. For large K, however, the convergence speeds of RCB with flat

ellipsoidal constraint and RCB with spherical constraint are about the same as that

of SCB; after convergence, the most accurate power estimate is provided by RCB

with flat ellipsoidal constraint.






























| -60 -60- | ;

-80 -80 -
-40 I I I -"

90 -60 -30 0 30 60 90 -190 -60 -30 0 30 60 90
0 degree 0 degree


- RCB (Flat ellipsoid) -- RCB (Flat ellipsoid)
RCB (Sphere) RCB (Sphere)
20 SCB 20 - SCB
20 -

a 0 .-40. '






-80o I 80 -

90 -60 -30 0 30 60 90 90 -60 -30 0 30 60 90
0 degree 0 degree

(c) (d)

Figure 2.8: Comparison of the beampatterns of SCB, RCB (with flat ellipsoidal
constraint) and RCB (with spherical constraint) when 6 = 2.40. (a) K = 1 and using
R. (b) K = 1 and using R with N = 10. (c) K = 8 and using R. (d) K = 8 and
using fR with N = 10. The true SOI power is 20 dB and eo = 3.1349 (corresponding
to A = 2).








38
























25 T 1 1 T r r 1 1 25 I 1 i I 1 1 T 1
2 .- --0-- -.-- ------o----o----

a 15 15

S10 10-
,,, l2 5- .. .-------'---'--'---'--"



-5 /

-10 -A RCB (Flat ellipsoid) -10o- RCB (Flat ellipsoid)
-- RCB (Sphere) -w RCB (Sphere)
S-o- SCB -a- SCB
10 20 30 40 50 60 70 80 90 100 10 20 30 40 50 60 70 80 90 100
Number of Snapshots Number of Snapshots

(a) (b)

Figure 2.9: Comparison of the SOI power estimates, versus N, obtained using SCB,
RCB (with flat ellipsoidal constraint) and RCB (with spherical constraint), all with
R when 6 = 2.4. (a) K = 1. (b) K = 8. The true SOI power is 20 dB and
co = 3.1349 (corresponding to A = 20).















CHAPTER 3
APPLICATION OF RCB TO AEROACOUSTICS
3.1 Introduction

The main motivation for our work in this chapter comes from an acoustic

imaging application in which the goal is to consistently estimate the SOI in the

presence of strong interference as well as some uncertainty in the SOI direction of

arrival. Due to its sensitivity to steering vector mismatch and small sample size, SCB

has not been used very much in acoustic imaging despite its potential benefits. The

various advantages of RCB, including robustness against array steering vector errors

and small sample size, high resolution, and superb interference suppression capability

make it a very promising approach to cure the problem of SCB.

Although RCB was devised under the narrowband assumption, it can also

deal with wideband acoustic signals by first dividing the array outputs into many

narrowband frequency bins using the fast Fourier transform (FFT) and then apply-

ing the narrowband RCB to each bin separately. However, it is well-known that as

the frequency increases, the beamwidth of both data-independent and data-adaptive

beamformers decreases. This beamwidth variation as a function of frequency will

subject the signals incident on the outer portions of the main beam to lowpass filter-

ing and lead to distorted signal spectra or inaccurate SOI power estimation [74]-[76].

Hence it is desirable that the beamwidth of a beamformer remains approximately

constant over all frequency bins of interest. In fact, a constant-beamwidth beam-

former is desirable in many applications including ultrasonics, underwater acoustics,

acoustic imaging and communications, and speech acquisition [74], [76], [77]. This

prevents future corrections for different frequencies, and contributes to consistent









sound pressure level (SPL) estimation, which means that for an acoustic wideband

monopole source with a flat spectrum the acoustic image for each frequency bin stays

the same.

In the past five decades, many approaches have been proposed to obtain

constant-beamwidth beamformers including harmonic nesting [75], [77]-[79], multi-

beam [74], [80], [81], asymptotic theory based methods [82], and approximation of a

continuously distributed sensor via a finite set of discrete sensors [83]-[85]. Among

these approaches, harmonic nesting is commonly used for acoustic imaging via mi-

crophone arrays. For example, in [77]-[79] a shading scheme is used for a directional

array consisting of a set of harmonically nested subarrays, each of which is designed

for a particular frequency bin. For each array element, shading weights are devised as

a function of the frequency. This shading scheme can provide a constant beamwidth

for frequencies between 10 and 40 kHz when used with the delay-and-sum (DAS)

beamformer. Hereafter, this approach will be referred to as the shaded DAS (SDAS).

In this chapter, we show that we can achieve a constant beamwidth across the

frequency bins for an adaptive beamformer, by combining our RCB with the shading

scheme devised in [77]-[79] provided that there are no strong interfering signals near

the main beam of the array. We refer to this approach as the constant-beamwidth

RCB (CBRCB) algorithm [86]. We also show that we can attain a constant pow-

erwidth, and hence consistent power estimates across the frequency bins, by using

RCB with a frequency-dependent uncertainty parameter for the array steering vector;

we refer to the so-obtained beamformer as the constant-powerwidth RCB (CPRCB)

[86]. CBRCB and CPRCB inherits the strength of RCB in the robustness against

array steering vector errors and finite sample size problems, high resolution, and

excellent interference suppression capability. Moreover, they both can be efficiently

implemented at a comparable computational cost with that of SCB.









3.2 Data Model and Problem Formulation of Acoustic Imaging

We focus herein on forming acoustic images using a microphone array, which

are obtained by determining the sound pressure estimates corresponding to the two-

or three- dimensional coordinates of a grid of locations. The signal at each grid

location of interest is referred to as the SOI.

First we introduce a wideband data model. Assume that a wideband SOI

impinges on an array with M elements. We divide each sensor output into N non-

overlapping blocks with each block consisting of I samples. We then apply an I-point

FFT to each block to obtain I narrowband frequency bins. The data vector, yi(n),

for the ith frequency bin and the nth snapshot can be written as

yi(n) =ai(xo)si(n)+ei(n), n = l, .,N; i = l, -,I (3.1)

where si(n) stands for the complex-valued waveform of the SOI that is present at the

location coordinate xo and the ith frequency bin, ei(n) represents a complex noise

plus interference vector for the coordinate xo and the ith frequency bin, and ai(xo) is

the SOI array steering vector, which depends on both xo and the ith frequency bin.

The nominal or assumed array steering vector ai(xo) has the form:
(ro [ Lej27rfir1 2r 1 eL i j27rfiTM (3.2)
6ij2nfi)- o r r2 rM (3

where rm and TO denote the propagation time from the source at xo to the mth sensor

and the array center, respectively. rm and ro represent the distance from the steered

location to the mth sensor and the array center, respectively. fi is the center of the

ith frequency bin. Note that this steering vector corresponds to the measurement at

the array center for the signal at the steered location. Due to the spherical spreading

of acoustic waves, the actual sound pressure at the steered location should be the one

at the array center multiplied by ro. If the steered location is in the far-field of the

array, we can use


a,(xo) = [ e2fi eij27fi2 ... ej27fiTM ]T


(3.3)









The covariance matrix for the ith frequency bin can be written as

Ri = E[yi(n)y;(n)], i= 1,...,I (3.4)

where E [.] is the expectation operator, and (.)* denotes the conjugate transpose.

In aeroacoustic measurements using arrays, the sound pressure response is

normally shown. The intensity of the sound pressure response is measured on a

logarithmic scale referred to as the sound pressure level (SPL), which is defined as

[87]
SPL = 20log1o(prms/prer) dB, (3.5)

where prm denotes the root-mean-squared pressure in Pa and Pref stands for the refer-

ence pressure. For air, the reference pressure is 20 pPa corresponding to the hearing
threshold of young people.

Next, we determine a scaling coefficient needed for calculating SPL estimates

based on R-. Let y(l), I = 1, 2,.. I, represent a wideband time-domain sequence
with I samples, and let Y(i), i = 1, 2, .. I, denote the frequency-domain sequence

obtained by applying an I-point FFT to {y(1)}. Let yi(l),l = 1,2,...,I, be the

narrowband time-domain sequence corresponding to the ith frequency bin of the
above frequency-domain sequence, in other words, {yi(l)} is obtained by using an
inverse FFT on {0, - 0, Y(i), 0, - 0}. Making use of only Y(i), we can determine
a root-mean-squared pressure estimate for the ith frequency bin, rm,, as follows:


Prm = i(1)l2 = jY(i)12 (3.6)
1=1

The second equality is due to the fact that =1 lyi(1)12 = lY(i)2 according to the
Parseval's theorem. Substituting (3.6) into (3.5), the SPL estimate for the ith nar-

rowband frequency bin can be written as 10 logio(JY(i)J2/(p22I2)) in dB. Therefore,

a scaling coefficient, -, should be used when estimating the sound pressure from

the data in the frequency domain. Let ai denote the root-mean-squared pressure









estimate for the ith narrowband frequency bin obtained using LRi. Then, according
to the previous discussion,

SPL = 10 loglo(oa/(p2reI2)) dB, (3.7)

In this chapter we distinguish between beampattern and powerpattern and
make use of this distinction in the design of the robust constant-beamwidth and
constant-powerwidth beamformers. We define the beampattern for the ith frequency
bin as

BPi(x) = la (x)w(x,)|2 (3.8)

where w(xg) denotes the beamformer's weight vector corresponding to a given loca-

tion, xg, and x is varied to cover each location of interest.
Next, we introduce the powerpattern, which at the ith frequency bin is defined
as
PPi(x) = |a!(x,)w(x)|2 (3.9)

We remark that the beampattern shows how the beamformer will pass the SOI and

interfering signals when it is steered to xg, whereas the powerpattern shows how the

beamformer will pass the signal at Xg when it is steered to x. PPi(x) can be used to
measure approximately the normalized power responses corresponding to a series of
locations, and hence it is named powerpattern.
To see this, we assume that the theoretical covariance matrix for the ith

frequency bin has the following form:

R, = a2ai(xo)a *(xo) + Q, (3.10)

where ao denotes the SOI power and Q stands for the interference-plus-noise covari-
ance matrix. If xg = xo and the signal-to-interference-plus-noise-ratio (SINR) is high,

then it follows that w*(x)Riw(x) aojlw*(x)ai(xg)|2 Oc PP,(x).
We next use an imaging example to illustrate the concept of beampattern
and powerpattern. Consider a function |a*(xl)w(z2) of two coordinate variables, xl









and x2. Then the beampattern is a slice of this function for x2 fixed, whereas the

powerpattern is a slice for xl fixed. Note that the beampattern and powerpattern of

the DAS beamformer are identical since its weight vector and steering vector have

the same functional form. However, this is not the case for an adaptive beamformer

due to the fact that its weight vector depends not only on the corresponding steering

vector, but also on the data.

Without loss of generality, we consider herein two-dimensional (2-D) array

imaging, in which the beamwidth or powerwidt is defined as the diameter of a circle

having the same area as the 3-dB contour of the main lobe of a 2-D beampattern or

powerpattern. The beamwidth shows how the nearby signals impact the estimation

of SOI. The powerwidth, on the other hand, shows how SOI impacts the estimation

of the nearby signals. From now on, we will concentrate on the ith frequency bin.

3.3 Constant-Powerwidth RCB

The beamwidth of RCB decreases with the frequency but generally it does not

depend on the choice of e if the SOI is far from the interference. On the other hand,

the powerwidth of RCB depends on the signal-to-noise-ratio (SNR), the frequency,

and e. Since the steering vector and hence its uncertainty set are both functions of

the frequency (see (3.2) and (2.10)), it is natural to consider altering the uncertainty

parameter e with the frequency. Intuitively, a larger e will yield a larger powerwidth.

By choosing a frequency-dependent parameter e for RCB, we obtain the constant-

powerwidth robust Capon beamformer (CPRCB), which is able to provide consistent

SPL estimates across the frequency bins for the source of interest [86]. However, the

beamwidth of CPRCB changes with the frequency in the same way as that of RCB.

Although it is difficult to yield an analytical formula for choosing e as a func-

tion of frequency to guarantee a nearly constant powerwidth, such a choice can be

readily made numerically via a contour plot of the powerwidths of RCB with respect









to the frequency and e. Given a desired powerwidth, we can determine e as a function

of frequency from the contour plot.

3.4 Constant-Beamwidth RCB

In [77]-[79] a shading scheme is used for a directional array consisting of a

set of harmonically nested subarrays, each of which is designed for a particular fre-

quency bin. For each array element, shading weights are devised as a function of

the frequency. This shading scheme can provide a constant beamwidth for frequen-

cies between 10 and 40 kHz when used with the DAS beamformer. We refer to this

approach as the shaded DAS (SDAS). Our RCB can be readily combined with the

shading scheme of [77], [78] to obtain a constant beamwidth for the desired frequency

band, as explained below. We refer to this approach as the constant-beamwidth RCB

(CBRCB) algorithm [86].

Let v denote the M x 1 vector containing the array element shading weights

for the given frequency bin. The assumed array steering vector for CBRCB can now

be written as:

ai = v i, (3.11)

where 0 denotes elementwise multiplication [5]. Accordingly, the covariance matrix

for the CBRCB is tapered as follows:

R 1 = Ri O (vvT). (3.12)

Since both R and vvT are positive semi-definite matrices, R is also positive semi-

definite [5], [88]. Note that RCB can be viewed as a special case of CBRCB with all

the elements of v being 1.

Similarly to (2.20) in Section 2.3.1, the CBRCB has the form:

max 02 subject to Ri a2aia* > 0
a2 -ai
Ilai i 112 < E, (3.13)









Table 3.1: Comparison of the features of CBRCB and CPRCB
Approach Constant Beamwidth Constant Powerwidth Loss of DOF
CBRCB Yes Yes Yes
CPRCB No Yes No


which can be solved like the RCB problem.

The ability of CBRCB to retain a constant beamwidth across the frequen-

cies makes it suitable for many applications such as speech acquisition. Furthermore,

CBRCB can also achieve constant powerwidth, which is essential for consistent imag-

ing. However, CBRCB is less powerful and flexible than CPRCB for applications

where constant powerwidth is demanded. First, CPRCB has more degrees of freedom

for interference suppression than CBRCB since at each frequency the shading scheme

involved in the latter tends to deactivate some elements of the full array. In addition,

CPRCB can be used with arbitrary arrays, while a special underlying array structure

is required for CBRCB due to the particular shading scheme employed. As men-

tioned earlier, SDAS can also be utilized to yield constant beamwidth. Nevertheless,

SDAS is data-independent and hence has poorer resolution and worse interference

suppression capability than CBRCB and CPRCB. The features of the CBRCB and

CPRCB approaches are listed in Table 3.1 in terms of constant beamwidth, constant

powerwidth and loss of degree of freedom (DOF). The computational costs of the

DAS, RCB, SCB, CBRCB and CPRCB approaches are listed in Table 3.2 in flops.

3.5 Numerical Examples

We provide some simulated examples to compare the performances of the

DAS, SDAS, SCB, RCB, CPRCB and CBRCB approaches for acoustic imaging. We

use the Small Aperture Directional Array (SADA) [77], [78], which consists of 33

microphones arranged in four circles of eight microphones each and one microphone

at the array center. The diameter of each circle is twice that of the closest circle










47











Table 3.2: Comparison of the computational cost of DAS, RCB, SCB, CBRCB and

CPRCB

Approach DAS SCB RCB CBRCB CPRCB

Computation O(M2) O(M3) O(M3) O(M3) O(M3)


0


3-
0 0
2- o
0 0
o o

000

-1C
%. 0 0 0 0 0 0 0 0

000
0 0
0
0 0
0
-2 -o

-3-
0
-4 -2 0 2 4
x (In)

(a)
4A. -,,


11 o
o o
oo
01 ooooo
0o0
o o
0 0
-1 -

-2-

-3-


-4 -2 0 2 4
x (In)

(b)

0
o
3-
0 0


o o
1 -

0-0 0 0 o 0

-1
0 0
-2 -

0 0
-3-
0


-4.
-4 -2 0 2 4 -4 -2 0 2 4
x (In) x (In)

(c) (d)


Figure 3.1: Microphone layout of the Small Aperture Directional Array (SADA) and

its three clusters. (a) SADA. (b) Cluster 1. (c) Cluster 2. (d) Cluster 3.


0
0 0
a 0
o o

0 0 0 0 0
0 0 0
0 o
0 0
0


A,




















f(Hz) x 10

Figure 3.2: Cluster shading weights for SADA, as functions of the frequency, with
wl, w2 and w3 corresponding to Cluster 1, Cluster 2 and Cluster 3, respectively.

it encloses. The maximum radius of the array is 3.89 inches. Figure 3.1 shows the

microphone layout of the SADA and its three subarrays used in the shading scheme

(referred to as clusters in [77], [78]). Note that some microphones are shared by

different clusters. Each cluster of SADA has the same directional characteristics for

a given wavenumber-length product kD,, where k is the wavenumber and D, is the

diagonal distance between the elements of the nth cluster. The wavenumber-length

products at 10 kHz for Cluster 3, at 20 kHz for Cluster 2, and at 40 kHz for Cluster

1 are the same. According to the array coordinate frame, the array is located in the

x-y plane, with center location at (0, 0, 0). Note that inch is used as the unit for the

3-D coordinates.

We assume that the distance between the array and the source is known and

plot the 2-D images by scanning the locations on a plane parallel to the array and

situated 5 feet above. We assume that a belongs to the uncertainty set


Ila- a112 (3.14)

where e is a user parameter chosen to account for the steering vector uncertainty.

Note that this form of uncertainty set used in the CPRCB can cover many kinds of

array errors, including calibration errors, look direction errors, or array covariance









estimation errors due to a small snapshot number (sample size). The uncertainty set

for the CBRCB is the same as above except that 5 is used instead of a. Figure 3.2

shows the SADA cluster shading weights as a function of the frequency bins, with wl,

w2 and w3 corresponding to Cluster 1, Cluster 2 and Cluster 3, respectively. Since

some array elements are shared by different clusters, the shading weights of those

elements are the sum of the corresponding cluster shading weights.

In the simulated examples below, we consider an array that is identical to

SADA, with a wideband monopole source (flat spectrum from 0 Hz to 70 kHz) located

at (0, 0, 60) (except for Figure 3.9) in the array coordinate frame and a spatially

white Gaussian noise with SNR equal to 20 dB and the SPL equal to 20 dB for each

frequency. We use an 8192-point FFT on the non-overlapping blocks (each containing

8192 samples) of simulated data to convert the wideband signal into 8192 narrowband

frequency bins.

In practice, the length for FFT operations is often chosen to be the nearest

power of two with respect to the number of samples per block. The minimum number

of frequency bins, which is required to decompose a wideband signal into narrowband

bins, can be determined by the following constraint [8]:

7ax < 1 (3.15)

where Tmax denotes the maximum time delay across the array aperture, B is the

bandwidth of each FFT bin, and < means much smaller than. It follows that

A 1
-- < (3.16)
c fs/I

where A denotes the array aperture, c is the sound speed, f, is the sampling frequency

and I represents the number of FFT bins. Inserting the array parameters, we obtain

I > 81. Hence, 8192 samples per block satisfies the above constraint. A smaller

sample number can also be used, e.g., 1024 samples per block with similar imaging

results.







50

20 M c
-0- CBRCB
18 -a- RCB
W- SDAS
16 DAS


12 -
6-




f (Hz) X.10
C I I f I I


Figure 3.3: Comparison of the beamwidths for the CBRCB, RCB, SDAS and DAS
methods with N = 64. We used e = 2.0 for RCB and CBRCB.


Example 3.1: Comparison of CBRCB, RCB, SDAS and DAS in

terms of beamwidth, powerwidth, and consistency of acoustic imaging.

Figure 3.3 compares the 3-dB beamwidths as functions of the frequency, correspond-

ing to the CBRCB, RCB, SDAS and DAS methods when N = 64. Note that herein

the beamwidth of RCB coincides with that of DAS and the beamwidth of CBRCB

coincides with that of SDAS. CBRCB and SDAS achieve constant beamwidth over

the frequency band from 10 to 40 kHz, whereas the beamwidths of RCB and DAS are

frequency dependent and decrease appreciably with the frequency. We used e = 2.0

for RCB and CBRCB. Other choices of e for RCB and CBRCB yield the same results

and hence they are not shown here. We remark that one can not achieve constant

beamwidth for RCB by varying e. For example, it can be shown that the beampat-

tern of RCB is independent of c if Q in (3.10) is white Gaussian noise, i.e., Q = oI,

where a, denotes the noise power and I is an identity matrix.

Figure 3.4 compares the 3-dB powerwidths of the CBRCB, RCB, SDAS and

DAS methods, as functions of the frequency, when N = 64. Note that the power-

widths of the DAS and RCB methods decrease drastically as the frequency increases,

while SDAS and CBRCB can both achieve approximately constant powerwidth. In

addition, CBRCB has much smaller powerwidth than SDAS. As can be seen from the







51


20C 2
-e- CBRCB -0- CBRCB
18 -- RCB 18 RCB
*- SDAS *a- SDAS
16- DAS 16- DAS
14 14 -


8-8
?12-\ ?12 .-
"). -. 0)



-..-... . ._ ._ . .
4- 4-L
2 1
C I I "_B_.__-B .-_--_-, -.__._._.__ 0- II_-- --.I
1 1.5 2 2.5 3 3.5 4 1 1.5 2 2.5 3 3.5 4
f(Hz) xo1 f (Hz) x 0o
(a) (b)

Figure 3.4: Comparison of the powerwidths for the CBRCB, RCB, SDAS and DAS
methods with N = 64 when (a) e = 1.0 for RCB and e = 0.5 for CBRCB and (b)
e = 2.0 for RCB and e = 1.0 for CBRCB.


figure, we can also adjust the powerwidth for CBRCB and RCB by choosing different

values of e.

Figure 3.5 compares the acoustic imaging results or sound pressure level (SPL)

estimates obtained via the DAS, SDAS, RCB and CBRCB methods for the narrow-

band frequency bins at 10 kHz and 40 kHz, with N = 64. The z axes show the

SPL. We used e = 2.0 for RCB and e = 1.0 for CBRCB. Note that we choose e for

CBRCB to be one half of that for RCB due to the fact that the squared norm of the

steering vector for CBRCB is about one half of that of RCB. As can be seen, the

DAS method has poor resolution and high sidelobes and its images vary considerably

with the frequency. RCB cannot be used to obtain consistent imaging results over

different frequency bins, either, though it has much better resolution than DAS. It

is worth noting that both SDAS and CBRCB maintain approximately the same SPL

estimates across the frequency bins, but the latter has much better resolution and

lower sidelobes and hence better interference rejection capability than the former.

It is obvious that CBRCB significantly outperforms the other methods. According

to the previous discussions and the results shown in Figures 3.3 and 3.4, it is the















20, -20,

2 1 2"2 2

(ft) -2 -2 x (t) -2 -2 X (ft)

(a) (b)

30 30, .






2 2 0 2
10 10,






2 x () x(t)
(c) (d)

Figure 3.5: Comparison of the acoustic imaging results obtained via the DAS, SDAS,
RCB and CBRCB methods with N = 64. (a) DAS with f = 10 kHz. (b) DAS with
f = 40 kHz. (c) SDAS with f = 10 kHz. (d) SDAS with f = 40 kHz. (e) RCB with
f = 10 kHz. (f) RCB with f = 40 kHz. (g) CBRCB with f = 10 kHz. (h) CBRCB
with f = 40 kHz. For RCB, e = 2.0. For CBRCB, E = 1.0. The z axes show the
SPL.

constant powerwidth rather than the constant beamwidth that contributes to the

better performance of CBRCB as compared to SDAS.

Example 3.2: Comparison of CPRCB and CBRCB in terms of con-

sistency of acoustic imaging. Figure 3.6 shows the contours of the 3-dB pow-

erwidth of RCB as e and the frequency vary, when N = 64. As can be seen, the

contours are almost linear with respect to the frequency and e. Therefore, given

a desired powerwidth, we can readily determine e as a function of the frequency

from the corresponding contour plot. Then CPRCB will have a constant powerwidth

across the frequency bins.


























































-1

) -2 -2


(g)


30,
20

10,

0
-10,
-20


0

y (f) -2


-10 1






o
-201
-30>
2
2



0()






030
30-,

20,

10 ..

m 0

-10-



-20

-30
21

0

x (ft)




Figure 3.5: Continued.


f(Hz) x10W


Figure 3.6: Contour plots of the powerwidth, versus e and the frequency, for the RCB

method. The numbers on the contours are the 3-dB powerwidths in inch.


2

01

-2 -2 x(t)


(h)


-2 x(


30,

20,

10-

o -^
-10

-20

-30>
2


-t
x











30 .
20-
10.
S00 ,
-10 1
-201


y (t)


(a) (b)

Figure 3.7: Acoustic imaging results obtained via the CPRCB method with N = 64.
(a) f = 10 kHz and e = 1.3. (b) f = 40 kHz and e = 13.



30 30 ..




-20, -20,



y( -1 -1-30
S -2 -2( -2 -2x()

(a) (b)

Figure 3.8: Acoustic imaging results obtained via the CBRCB method with N = 64
and e = 0.65. (a) f = 10 kHz. (b) f = 40 kHz.


In Figure 3.7, we show the imaging results obtained via the CPRCB approach

by choosing e = 1.3 when f = 10 kHz and e = 13 when f = 40 kHz from the 3 inch

powerwidth contour in Figure 3.6. The similarity of the SOI SPL estimates obtained

with CPRCB at these two frequencies, especially near the powerwidth area, verifies

the consistency of CPRCB in powerpattern across the frequencies.

Figure 3.8 shows the imaging results obtained via the CBRCB approach with

e = 0.65, for f = 10 kHz and f = 40 kHz. Again we note the consistency in the

imaging results. Therefore, both CPRCB and CBRCB are suitable for applications

where consistent SPL estimates are desirable. However, the sidelobes in Figure 3.7(b)










I BIR-
-U- CPRCB -0- CPRCB
40- SCB 35 -0- SCB .
-A- SDAS -A- SDAS ..- *
D- AS -- DAS
30 30-
O 20 -n -..,- ,,.A 5- --Ai








(a) (b)

Figure 3.9: Comparison of the SINR and SOI SPL estimates obtained via the
CBRCB, CPRCB, SCB, SDAS and DAS methods, versus the number of interfer-
estimate. For CPRCB, = 2.0. For CBRCB, 1.0. We consider a look direction






error case where the assumed source location is (0, 0, 60) but the actual point source
is located at (0.2, 0.2, 60) with SNR equal to 20 dB. The INRs are equal to 40 dB.
K K






are higher and rougher than those in Figure 3.8(b). Despite this fact, CPRCB does

Figunot perform worse thamparison CBRCB, see SINR and SOI SPL estimates obtained via the
CBCB,Example 3.3: Comparison of CPRCB, CBRCB, SCB, SDAS and DAS methods, versus the number of interfer-
ences K, for the narrowband frequency bin at f = 20 kHz. (a) SINR. (b) SOI SPL






DAS in the presence = 2.0.of look direction e = errors. We consider a look direction
error case where the assumed source location is (0, 0, 60) but the actual point source
located at (0.2, 0.2, 60) with SNR equal to 20 dB. Also we consider are equal to 40 dB.
are higher and rougher than those in Figure 3.8(b). Despite this fact, CPRCB does

not perform worse than CBRCB, see the next example.

Example 3.3: Comparison of CPRCB, CBRCB, SCB, SDAS and

DAS in the presence of look direction errors. We consider a look direction

error case where the assumed source location is (0, 0, 60) but the actual source is

located at (0.2,0.2,60) with SNR equal to 20 dB. Also we consider a varying number

of interference from K = 0 to K = 20, which are situated on a circle with a radius

of 20 inches and have an INR equal to 40 dB. The circle is on a plane parallel to

the array and situated 60 inches above. We assume that the theoretical covariance

matrix R is known here.

Figure 3.9 compares the SINR and SOI SPL estimates obtained via the CBRCB,

CPRCB, SCB, SDAS and DAS methods, versus the number of interference K, for

the narrowband frequency bin at f = 20 kHz. For CPRCB, e = 2.0. For CBRCB,

e = 1.0. Note that SCB is very sensitive to the steering vector mismatch and suffers

from severe performance degradation in SINR and SOI SPL estimates. Although DAS


501 1 I


4U









and SDAS are robust against array errors, they have poor capacity for interference

suppression. Consequently, their SINRs and SOI SPL estimates are unsatisfactory.

CBRCB and CPRCB outperform the other approaches due to their robustness to

steering vector errors, better resolution and much better interference rejection capa-

bility than the data-independent beamformers. As can be seen, CPRCB has higher

SINR than CBRCB. This is due to the fact that the former has more degrees of free-

dom (DOFs) for interference suppression than the latter. It might seem surprising

that the performance of SCB improves as the number of interference K increases.

There is a simple explanation for this. When K is small, SCB has enough many

DOFs and the SOI is suppressed as interference. As K increases, SCB focuses more

on suppressing the interference than the SOI since the INR is much higher than the

SNR.















CHAPTER 4
APPLICATION OF RCB TO ULTRASONICS
4.1 Introduction

The data-independent delay-and-sum (DAS) beamformer is widely used for

ultrasound imaging applications [89]-[91], although it has lower resolution and worse

interference suppression capability than an adaptive beamformer such as the standard

Capon beamformer (SCB) [10], provided that the array steering vector corresponding

to the signal of interest (SOI) is accurately known. However, there are many fac-

tors in practice that can degrade the performance of SCB including pointing errors

due to inaccurate knowledge of the source location, array calibration errors due to

imprecise knowledge of the transducer responses and positions, and array covariance

matrix estimation errors due to a small sample size. Owing to these problems, the

performance of SCB may become worse than that of the standard data-independent

beamformers such as the DAS beamformer. Consequently, SCB has not been used

extensively despite its potential benefits.

Much work has been done to improve the robustness of SCB over the past

three decades and the literature on robust adaptive beamforming is extensive. How-

ever, most of the early suggested methods do not directly address the problem of the

uncertainty of the array steering vector [51], [53]-[55], [58], [60]. The robust Capon

beamformer (RCB) we presented in [13], [14], on the other hand, couples the formu-

lation of the covariance fitting based SCB in [12] with an ellipsoidal uncertainty set

of the array steering vector. In addition, RCB comprises a simple way of eliminating

the scaling ambiguity when estimating the power of the desired signal.









The aforementioned approaches are all devised based on the narrowband as-

sumptions. Wideband signals are often encountered in ultrasound imaging applica-

tions such as pulse-echo detection. For a wideband signal, we can divide the array

outputs into many narrowband frequency bins using the Fourier transform, apply the

narrowband beamformer for each frequency bin and then combine the narrowband

estimation results. However, such incoherent approaches have the disadvantage that

their performances are limited by the threshold effects of the observation time and the

signal-to-noise ratios of the individual narrowband beamformers [92], [93]. To make

better use of the large bandwidth, several methods have been proposed for wideband

far-field processing using passive arrays. For example, the coherent signal-subspace

method (CSM) in [94], [95] utilizes focusing matrices to combine the signal subspaces

at different frequencies into one focused covariance matrix. However, the estimation

errors of the focusing matrices may result in the source location bias. Such bias can

be avoided by the steered minimum variance (STMV) method in [96], which uses a

steered covariance matrix in the space-time domain, obtained after inserting proper

time delays to the data samples for each direction of arrival. Nevertheless, neither

CSM nor STMV can be directly applied to ultrasound imaging, which often requires

wideband near-field processing using active arrays.

To deal with the ultrasound imaging problem, we may explore two options.

The first one borrows the idea from the STMV method in [96] and the observation that

the number of transmitters in active arrays corresponds to the number of snapshots

in passive arrays [97]. By aligning the received signals from a complete multistatic

data set to the focal point based on the appropriate time delays, we can construct a

pseudo-covariance matrix for beamforming. The second one comes from exploiting a

covariance matrix interpretation of the time reversal operator as shown in [97]. In the

past decade, time reversal has been a very active domain of research and has found

many applications including medical diagnostics and therapeutics, nondestructive









evaluation, underwater acoustics, acoustic room dereverberation, communications

[98]-[104]. Time reversal can be cast into two main categories as physical time reversal
and synthetic time reversal. For the physical time reversal, the recorded signals at

the array are time reversed and physically re-emitted into the medium to optimally

focus on the targets, whereas for the synthetic time reversal, the backpropagation is

done numerically to form images. The essential distinction between the physical time

reversal and synthetic time reversal is that the former does not require knowledge of

the medium, while the latter does. Specifically, the Green's function in the medium is

needed for synthetic time reversal. In practice, the Green's function can be measured

as used for acoustic room dereverberation [99], [105]. Also distorted Born iterative

(DBI) method and the multiple frequency DBI method can be used to estimate the

inhomogeneous Green's function in the medium [106], [107].

In this chapter, we show that RCB can be used for ultrasound imaging in

two ways, one based on time delay and the other based on time reversal. The time-

delay based RCB utilizes a pseudo-covariance matrix constructed in the space-time

domain similarly to [96], while the time-reversal based RCB exploits a covariance

matrix interpretation of the time reversal operator in the frequency domain as shown
in [97]. Due to its inherent robustness, RCB can allow array calibration errors,

caused by transducer magnitude and phase response mismatch as well as transducer

position errors, and array covariance matrix estimation errors, as a result of small

sample sizes. In addition, the time-delay based RCB can tolerate the misalignment

of data samples and the time-reversal based RCB can withstand the uncertainty of

the Green's function. Furthermore, RCB has the desirable features including high

resolution and low sidelobes, both of which are very important for ultrasound imaging

[91], [108], [109]. Due to its adaptivity, RCB has higher resolution than DAS. RCB
belongs to the class of the (extended) diagonal loading based approaches and it can

obtain the optimal diagonal loading level based on the uncertainty set of the array









steering vector [13]. Hence RCB has low sidelobes even in the presence of array errors

and small sample sizes, which characterizes the diagonal loading based approaches

[55]. RCB can also be efficiently computed at a comparable cost with that of SCB by

using the Lagrange multiplier methodology [13]. Finally, simulated and experimental

examples are presented to illustrate the effectiveness of RCB for ultrasound imaging.

4.2 Problem Formulation

Consider an active array of M transducers, which uses the pulse-echo mode to

probe the unknown propagating medium, i.e., pulses are emitted and data sequences

of the back-scattered echoes are stored. Depending on how data are acquired, there

are two major methods, namely a monostatic approach and a multistatic approach. In

a monostatic approach, a single array element is used and it acts both as a transmitter

and a receiver, whereas in a multistatic approach, a transmitter is scanned across the

array aperture with multiple receivers.

A complete data set, obtained by a multistatic approach for data acquisition,

is considered in the chapter. This means that the resulting multistatic array response

matrix P(t) with t denoting the time instant, includes the A-scan data for all possible

transmitter and receiver pairs of the array. Specifically, the A-scan data sequence

Pi,(t) is obtained by recording the back-scattered echo at the jth transducer, after
transmitting a pulse from the ith transducer. Assume that each A-scan data sequence

contains N data samples. Then P(t) has a total of M x M x N data samples.

Let x and z represent the lateral and axial axes of a coordinate system, respec-

tively, (xi, zi) denote the location of the ith transmitter, and (xj, zj) be the location

of the jth receiver. Let p = (XF, ZF) represent the location of a focal point of interest

in the imaging region. The time delay due to the ultrasonic wave propagation from

the ith transmitter to the focal point p and then back to the jth receiver is
1 1
Ti,j(P) = [(Xi XF)2 + ( Z F)2 [(Xj -F)2 ( ZF)2] (4.1)
c c









where c is the sound propagation speed in the medium.
Given a focal point p, we define yij(p) as the sample selected according to the
time delay Tj(p) from the A-scan data sequence Pi,j(t), i.e., yi,j(p) = Pi,j(Tj(p)).
Due to the acoustic impedance variations in the medium, either reflection or scatter-
ing occurs depending on the size of the target. We can write yi,j(p) as the sum of
two terms:


i,j(p) = A(p) + ei,j(p), i = M; j = 1, M (4.2)

where i3(p) denotes the on-axis signal that is proportional to the ultrasonic reflectivity
or scattering strength at the focal point p, and ei,j(p) denotes the residual term due to
the contributions from the noise and interference (off-axis signals) at other locations.
Let

Yi(p) = [ yi,I(P) Yi,2(P) Yi,M(P) ]T, i= 1,.. ,M (4.3)

and

ei(p) = [ ei,l(p) ei,2(p) '. ei,M(p) IT, i= 1,.-,M (4.4)

It follows that

Yi(p)=ai3(p)+ei(p), i=1, ..,M (4.5)

where yi(p) represents the aligned array data vector due to the ith transmitter, a
denotes the array steering vector and is approximately equal to 1MxI, which is an M
by 1 vector with unit elements. In practice the elements of the steering vector a are
also affected by data misalignment or array calibration errors due to the transducer
magnitude and phase response mismatch as well as transducer position errors.
Let Y(p) be a matrix with yi(p) as its ith column, i.e.,


Y(p) = [ yl(p) y2(P) '.. YM(P) ]


(4.6)









Then the pseudo-covariance matrix relating to the ultrasonic reflectivity or scattering

strength at the focal point p can be written as

R(p) = Y(p)Y*(p) (4.7)

where (.)* denotes the conjugate transpose. With probability 1, R(p) has full rank.

Note that R(p) is calculated in the same way as the sample covariance matrix for pas-

sive arrays. We refer to the resulting space-time matrix R(p) as a pseudo-covariance

matrix due to the fact that different columns of Y(p) correspond to different trans-

mitters and hence contain the contributions from the targets on distinct arcs, which

coincide on the same focal point.

Both SCB and the DAS methods can be applied to the data model in (4.5) for

ultrasound imaging. The goal of time-delay based SCB is to pass the on-axis signal

(the signal on the focal point) without distortion, and at the same time minimize the

contributions from the off-axis signals (other signals on the arcs). However, in practice

the time-delay based SCB suffers from severe performance degradations due to the

small sample size (note that the number of snapshots N is equal to the number of

transmitters M) and array steering vector errors caused by data sample misalignment

or array calibration errors. In fact, there is a close relationship between the cases

of steering vector errors and small-sample errors [38] in the sense that the difference

between the sample covariance matrix (estimated from a finite number of snapshots)

and the corresponding theoretical covariance matrix can be viewed as due to steering

vector errors. The conventional DAS method, on the other hand, simply adds up

all the signals whether they are on-axis or off-axis. It cannot adapt to the incoming

data and hence has poor resolution and off-axis signal suppression capability.

From the above discussions, we see that the time-delay based ultrasound imag-

ing in the presence of steering vector mismatch and small sample size needs to be

formulated as a robust adaptive beamforming problem. In particular, we will deal with

the following problem in this chapter: extend the standard Capon beamformer so as









to be able to determine accurately the ultrasonic reflectivity or scattering strength

for each focal point of interest, even when only an imprecise knowledge of its steering

vector is available. More specifically, we assume that the only knowledge we have

about a is that it belongs to the following uncertainty sphere:


{a I a all2 < e} (4.8)

where a and e are given (a is the assumed steering vector, and e is a user parameter

describing the uncertainty of the array steering vector). From now on, we will con-

centrate on the focal point p. For notational simplicity, whenever convenient, we will

drop the dependence on p in the following sections.

Next, we present some preliminaries on time reversal. From the array response

matrix P(t) in the time domain, we can readily obtain the array response matrix

P(w) in the frequency domain, which is symmetric but not Hermitian. Assume that

there are K isotropic point targets located at yi, ..., YK, without multiple scattering

among them and assume that the transducers are all isotropic point transmitters and

receivers. Then P(w) has the form [103]:
K
P(w) = f(w) ( (i(w)g(yi, w)g(yi, ) (4.9)
i=1

where f(w) is the Fourier transform of the transmitted pulse f(t), (i(w) denotes the

scattering coefficient of the ith target, (.)T denotes the transpose, and g(yi, w) is the

illuminating Green's vector onto the array from the point yi, which can be written

as [103]:


g(yi,w) = [ G(y,,xl,w) G(yi,x2,w) ... G(y, xM, w) ]T (4.10)

where xj denotes the location of the jth transducer, G(yi, xj, w) is the Green's func-

tion in the medium at the frequency w. For the homogeneous background in two-

dimensional free space with a known sound speed c, the Green's function can be









written as
-i
G(yi, x, w) = -4Ho (k||iy xjll) (4.11)

where k denotes the wavenumber (k = w/c), Ho is the zeroth order Hankel function

of the first kind. Note that in practice, P(w) in (4.9) should also contain an additive

term due to the noise.

The time reversal operator is defined as

T(w) = P*(w)P(w) (4.12)

As stated in [103], the time reversal operator contains information obtained from

probing the medium twice. The connection between the eigenvectors of the time

reversal operator and the scatterers was investigated and exploited by the D.O.R.T

method (French acronym for decomposition of the time reversal operator) [110]. Re-

cently it was further shown in [97] that the time reversal operator for active arrays

can be interpreted as a covariance matrix for passive arrays, and the illuminating

Green's vector in the medium can be used as the assumed steering vector. Therefore,

we can explore some algorithms originally developed for passive arrays for time-

reversal based imaging. The limitation, however, is that we need to have sufficient a

priori knowledge of the medium to obtain the Green's vector. If the knowledge is im-

precise, there will be steering vector mismatch and consequently, the performance of

a time-reversal based adaptive beamformer may degrade severely. This motivates us

to consider a robust adaptive beamformer. Similar to the time-delay based method,

we can apply RCB by allowing the Green's vector to be within an uncertainty set as

in (4.8).

4.3 Time-Delay Based RCB

In this section we focus on the time-delay based RCB, with additional dis-

cussions on time-delay based SCB and DAS. We show how to calculate the weight

vectors of these beamformers and use coherent processing for ultrasound imaging.









Note that for the time-delay based methods, we use R in (4.7) as the covariance

matrix and 1Mxl as the assumed steering vector a, i.e., a = 1Mxl.

The common formulation of the adaptive beamforming problem that leads to

the standard Capon beamformer, when the array steering vector is assumed to be a,

is as follows [10], [5].

(a) Determine the M x 1 vector w that is the solution to the following linearly

constrained quadratic problem:

minw*Rw subject to w*a = 1 (4.13)
w

(b) Use w*Rw as an estimate of the sound pressure u2.

The solution to (4.13) is easily derived:

R-1a R--1Mxl
S=- = (4.14)
a*R-la l1xlR-11Mxl

Using (4.14) in Step (b) above yields the following estimate of a2:

-2 = 1 1(4.15)
a*R-la 1MxlR-11Mx1

The problem with the standard Capon beamformer is that it is very sensitive

to steering vector errors and small sample sizes. In our recent papers [13], [14], we

have proposed a robust Capon beamformer (RCB), which extends the SCB so as to

be able to accurately determine the power of the signal-of-interest (SOI) even when

only an imprecise knowledge of its actual steering vector is available.

To derive the robust Capon beamforming approach, we use the reformulation

of the Capon beamforming problem in [12], to which we append the uncertainty set

in (4.8) [13], [14]. Proceeding in this way we directly obtain a robust estimate of a2,

without any intermediate calculation of a vector w:

maxo 2 subject to R a2aa* > 0
2 ,a
Ila- all12 < E (4.16)









Note that the first line above can be interpreted as a covariance fitting problem:

given R and a, we wish to determine the largest possible SOI term, a2aa*, that can

be a part of R under the natural constraint that the residual covariance matrix is

still positive semidefinite.

For any given a, the solution &2 to (4.16) is given by the counterpart of (4.15)

with a replaced by a, as shown in [14]. Hence (4.16) can be reduced to the following

problem

mina*R-la subject to Ila a2 < (4.17)
a
To exclude the trivial solution a = 0 to (4.17), we assume that

|II|a 2 > E (4.18)

Because the solution to (4.17) (under (4.18)) will evidently occur on the boundary of

the constraint set [13], we can re-formulate (4.17) as the following quadratic problem
with a quadratic equality constraint:

mina*R-la subject to la a112 = e (4.19)
a

The solution to the above optimization problem can be determined by using the

Lagrange multiplier methodology and it has the form [13], [14]:

a -= + I (4.20)

= (I + AR) (4.21)

where we have used the matrix inversion lemma [5] to obtain the second equality.

The Lagrange multiplier A > 0 can be obtained as the solution to the constraint

equation:

(I + AR)-1 l2 = 0 (4.22)

the left-hand side of which is a monotonically decreasing function of A for A > 0.

Hence A can be solved easily, e.g., by a Newton's method. We then use the so-
obtained A in (4.21) to obtain a. However, the norm of the resulting a tends to be








small due to the scaling ambiguity [13], [14]. Hence, we use the scaled estimate of

the array steering vector as follows:

a = ViMa/|lla (4.23)

The scaling above is based, without loss of generality, on the constant norm assump-

tion ||ai12 = M.

To obtain the weight vector for RCB, A in (4.23) is used to replace a in (4.14):

R-lA
w = (4.24)
a*R-1f

We remark that the major computational demand of our RCB comes from the eigen-

decomposition of the Hermitian matrix R, which requires 0(M3) flops (we can easily

inverse R in (4.24) based on the eigen-decomposition). Therefore, the computational

complexity of RCB is comparable to that of SCB.

When the pseudo-covariance matrix R in (4.13) is replaced with an identity

matrix corresponding to the pseudo-covariance matrix for spatially white noise, it

becomes a delay-and-sum beamformer as follows:

minw*w subject to w*a = 1 (4.25)
w

The solution is
w = 1M (4.26)
M M
From the weight vector in (4.26), we see that DAS is a data-independent approach.
Once the weight vector of a beamformer is obtained, we can estimate 3i as

(we reinstate the dependence on p for clarity):

i(p) = *(p)yi(p), i= 1, ,M (4.27)

We can then use coherent processing to compute the final estimate 3(p) as follows:
M
w(p) = M i(P) (4.28)
i=1









When {/(p)} are calculated for a grid of points in the imaging region of interest, we

obtain an ultrasound image.

Using (4.26) into (4.27) and then into (4.28), we obtain the following well-

known form of the DAS estimate based on a complete multistatic data set:
M M
() = y,(p) (4.29)
i=1 j=1

If instead a monostatic data set is used for DAS, it is easy to show that
M
M(p) = Y, (P) (4.30)
i=l

Note from (4.29) and (4.30) that the DAS using a multistatic data set exploits more

data samples than the DAS using a monostatic data set. Hence the image obtained

by the former should have better quality than that obtained by the latter at the

expense of extra data collection and more computational load [89]. This will also be

verified by our examples later on.

4.4 Time-Reversal Based RCB

In this section we discuss the time-reversal based RCB and several other al-

gorithms including the time-reversal based SCB, matched field (MF), and multiple

signal classification (MUSIC).

We consider herein a coherent approach, which exploits the time-bandwidth

product of the wideband signal. We use a scheme to construct the focusing matri-

ces for the time reversal operators, similar to the rotational signal subspace (RSS)

focusing approach in [95]. The difference is that the latter is used for passive arrays

in the far field and requires the initial estimates of the source locations, while the

former is used for active arrays in the near field and does not demand the initial

estimates of the source locations. Specifically, in the former method we use focusing

matrices to transform the time reversal operators at different frequency bins onto a









single reference frequency, and then sum the transformed time reversal operators to
obtain the coherently focused time reversal operator.
Let
Y(w) = P*(w) (4.31)

Substituting (4.31) into (4.12) yields

T(w) = Y(w)Y*(w) (4.32)

Assume that the reference frequency is equal to wo, and we wish to find the
unitary focusing matrices for J frequency bins at wj, j = 1, 2,..., J. Unitary focusing
matrices are desired due to the fact that they have better performances in preserving
the shape of the noise covariance matrix and the array SNR than nonunitary focusing

matrices. In [95], the following constrained minimization problem is used to determine
the unitary focusing matrices F(wj):

min |Y(wo) F(wj)Y(wj)IF subject to F*(wj)F(wj) = I (4.33)
F(wj)

where I|| IF denotes the Frobenius matrix norm. Let the singular value decomposition

(SVD) of Y(wj)Y*(wo) be:

Y(wj)Y*(wo) = U(wj)E(wj)V*(Wj) (4.34)

It can be proved that the solution to the problem in (4.33) is [95]

F(wj) = V(wj)U*(wj) (4.35)

Then the coherently focused time reversal operator is given by
J
T(wo) = F F(wj)T(wj)F*(wj) (4.36)
j=1

Note that for the coherent approach, T(wo) is used as the covariance matrix and the
illuminating Green's vector g(wo,p) is used as the assumed steering vector for the
focal point p.









Next, we briefly present several time-reversal based methods, whose formula-

tions are analogous to the time-delay based counterparts in the previous section. To

simplify the notation, we omit the argument wo in T(wo) and g(p, wo).

The time-reversal based RCB has the form:

mina*T-la subject to Ila g(p)|12 < e (4.37)
a

Note that the norm square of g(p) depends on the location of the focal point and is

not equal to M. Hence e should be chosen as a function of g(p). In addition, after

the estimated steering vector a is obtained, we now need to scale it based on the

constant norm assumption I|alI = llg(p)ll as follows:

a = II(p)lla/l|al| (4.38)

The remaining steps of the time-reversal based RCB are in direct analogy to those

of the time-delay based RCB, and are hence omitted herein for brevity.

The time-reversal based SCB can be written as:

minw*Tw subject to w*g(p) = 1 (4.39)
w

The solution is

w= (4.40)
g*(p)T-Ig(p)
The time-reversal based MF has the form:

minw*w subject to w*g(p) = 1 (4.41)
w

The solution is

W = () (4.42)
II (p ) I2I
As its name implies, the time-reversal based MF matches its weight vector to the

field of a virtual medium with the Green's vector equal to g(p) at the focal point p.

Note that the formulation of MF in (4.41) is similar to that of the time-delay based

DAS in (4.25).









Once we have obtained the weight vector w for the time-reversal based SCB,
RCB and MF, we can estimate the magnitude square of the on-axis signal as (we
reinstate the dependence on p for clarity):


I12(p) = w*(p)Tw(p) (4.43)

When {I/112(p)} are calculated for a grid of points in the imaging region of interest,
we obtain an ultrasound image.
In [102], MUSIC was considered for time-reversal based imaging. MUSIC
was a super-resolution approach originally proposed for passive arrays [17]. In the
following, we briefly review the MUSIC algorithm. We will compare its performance
with our proposed methods in Section 4.5. Let {vm}m= be the eigenvectors of T
associated with the eigenvalues arranged in a decreasing order, and let S and G
denote the following two matrices formed from {vm})M=:

S = [V1 * VK], G = [VK+ * VM] (4.44)

where K denotes the dimension of the signal subspace, e.g., the number of point
targets in the medium. Note that S and G satisfy

S*S= I, G*G = I and S*G = 0 (4.45)

The locations of the K targets are given by the positions of the K highest peaks of
the function
1
h(p) = *() (4.46)
g*(p)GG*g(p)
As is well-known, MUSIC is not an intensity estimator but a locator. It cannot cope
with distributed targets. Also MUSIC demands the information of the dimension of
the signal subspace K, which is unknown in practice and may be difficult to accurately
estimate in the presence of colored noise.









4.5 Simulated and Experimental Examples

In this section, we provide several simulated and experimental examples to
compare the performances of the time-reversal and time-delay based methods includ-

ing DAS, SCB, RCB, MF and MUSIC.

The Field II software [111] is used as a simulation tool to generate a complete
multistatic data set. In the simulations, the probing pulse is differential Gaussian
and has the form:


f(t) = -27r22 (t - exp [-2 2 (t ] (4.47)

where v denotes the central frequency (v = 3 MHz). The sampling frequency is 100
MHz. The array comprises 21 transducers with inter-sensor spacing (referred to as
pitch in ultrasound imaging) equal to A, where A is the central wavelength. The

transducers are isotropic and located at

xj = (-5A + A, 0), = 1,...,21. (4.48)

There are two point targets located at (-2A, 20A) and (3A, 19A) in a free-space homo-
geneous medium with a sound speed equal to 1540 m/s. In all the simulated images
below, a 30 dB display dynamic range is used with the plus symbols denoting the
true positions of the two targets. The unit for the axes is A. For the time-delay based
methods, each A-scan data sequence is preprocessed by a matched filter based on
the probing pulse. For the time-reversal based methods, we focus the time reversal
operators at all of the frequency bins onto the central frequency (corresponding to the
spectral peak of the transmitted pulse) to obtain a coherently focused time reversal
operator. Unless specifically clarified, the complete multistatic data set is employed.
Figures 4.1 and 4.2 show the imaging results obtained via the time-reversal and
time-delay based methods, respectively, for the case of no transducer position errors.

We see from Figure 4.1 that the time-reversal based MF has poor resolution and high









sidelobes. Although the time-reversal based SCB image has peaks corresponding to

the targets, it also contains some sidelobes like that of MF. Both MUSIC and the

time-reversal based RCB have good cross-range resolution and low sidelobes. On the

other hand, their range resolutions are not as good. This is due to the fact that the

Green's vector is more sensitive in the cross range than in the range. We can verify

this by imaging the norm of the difference of the Green's vector at a given point

with those at its surrounding points. To improve the range resolution of the time-

reversal based methods, one method is to combine the direction of arrival estimation

with the time delay estimation as in [103], [104] for imaging in the random media.

However, this method requires very complicated procedures to estimate the time

delay. Note from Figure 4.1 that the time-reversal based RCB outperforms the other

three methods. From Figure 4.2, we see that the time-delay based SCB yields a very

poor image. This is probably due to the high sidelobes caused by the small sample

size. Comparing the DAS image obtained using monostatic data and that obtained

using multistatic data, we observe that the targets are more localized in the latter

with weaker trailing effects near the targets. A simple explanation is that averaging

by using more samples lowers the noise contributions to the estimate at the focal

point and hence higher SNR will result. This is actually why we concentrate on the

multistatic data processing in this chapter. In terms of resolution and sidelobes, the

time-delay based RCB has the best imaging result among all the time-delay based

methods considered here. As can be seen from Figures 4.1 and 4.2, the time-delay

based approaches have better range resolution but worse cross-range resolution than

the time-reversal based counterparts.

Figures 4.3 and 4.4 are similar to Figures 4.1 and 4.2, except that now each

transducer position is perturbed by a zero-mean Gaussian random variable with vari-

ance equal to 0.01A/2. The perturbing Gaussian random variables are independent

of each other. This will lead to data misalignment for the time-delay based methods









and uncertainty in the Green's function for the time-reversal based methods. Hence

there are array steering vector errors. As a result, degradations occur in the ultra-

sound images. Note that the time-reversal based RCB and the time-delay based RCB

still achieve good performances, outperforming the other methods considered. This

verifies the robustness of RCB against steering vector errors.

Next, we present some experimental results to compare the performances of

the time-delay based DAS, SCB and RCB approaches. We did not consider the time-

reversal based methods here due to lack of information for the Green's function in

the experimental media. The complete multistatic data sets were obtained from the

Bioacoustics Research Laboratory of the University of Illinois at Urbana-Champaign

(the data set for the wire phantom in a complex pattern and the data set for the

rat mammary tumor) and the Biomedical Ultrasonics Laboratory of the University of

Michigan at Ann Arbor (all the other data sets). To reconstruct an ultrasound image,

we first use a digital bandpass filter on the data. Then we apply the aforementioned

beamforming algorithms on the filtered data. Finally, we apply envelope detection,

gain compensation, scan conversion and logarithmic compression. Except specifically

stated otherwise, all the images below are shown over a 50 dB display dynamic range.

We first consider a wire phantom containing six wire targets. The data are

collected using a 128-element linear array (M = 128). The transducer center fre-

quency is 3.5 MHz and the sampling rate is 13.8889 MHz. The sound velocity is

1480 m/s. In Figure 4.5, we show the ultrasound images obtained via DAS, SCB and

RCB using the entire array and only the central 64 elements of the array, respectively.

The latter can be regarded as using a smaller array with M = 64. Note that SCB

yields very poor images where the wire targets are lost. Although the wires are more

recognizable in the DAS images than in the SCB images, the sidelobes of the DAS

images are still very high. Clearly, RCB significantly outperforms DAS and SCB in

terms of resolution and sidelobes.









To take a closer look at the resolution of DAS and RCB, in Figure 4.6 we

compare the ultrasound images obtained via DAS and RCB using only the central

32 and 64 elements of the array over a 20 dB display dynamic range. This can be

regarded as using smaller arrays with 32 and 64 elements, respectively. As can be

seen from Figures 4.6(a) and 4.6(b), the resolution of the DAS improves if more

transducers are used. In Figures 4.6(c)-4.6(f), we vary e from 0.2 to 20 for RCB and

show the resulting images. Note that as e increases, the resolution of RCB becomes

lower. This is not surprising since e is a parameter used to adjust the robustness of

RCB. In fact, we can prove that the bigger the e, the smaller the Lagrange multiplier

A and the bigger the loading level. Consequently, the less significant role the pseudo-

covariance matrix R plays and the closer RCB is to DAS. On the other hand, the

smaller the e, the closer RCB is to SCB and the more sensitive RCB becomes to

the steering vector errors. From Figure 4.5, we can easily observe that, for a large

range of choices for e, RCB using an array with 32 transducers can achieve a similar

resolution as that of DAS using an array with 64 transducers, outperforming that

of DAS using an array with 32 transducers. This means that RCB can double the

resolution of DAS without complicating the array hardware at all.

As a further illustration, we use a data set containing some wire targets in

a more complicated pattern. The data are collected using a 64-element linear array

(M = 64). The transducer center frequency is 2.6 MHz and the sampling rate is

25 MHz. The sound velocity is 1450 m/s. In Figure 4.7, we show the ultrasound

images obtained via DAS, SCB and RCB. For DAS, we compare its images obtained

using monostatic data and multistatic data. Note that by choosing the A-scans with

each transducer acting both as a transmitter and the only corresponding receiver,

the monostatic data can be easily obtained from the multistatic data. As can be

seen from Figures 4.7(a) and 4.7(b), the image obtained based on multistatic data

is much better than that obtained based on monostatic data. Due to its sensitivity









to steering vector errors and small sample sizes, SCB yields a very poor image as in

Figure 4.7(c), in which the wires are hardly recognizable. The RCB image, on the

other hand, has very good resolution and low sidelobes, surpassing both SCB and

DAS in performance.

To investigate the contrast resolution, we consider now a cyst phantom con-

taining four cyst targets. The data are collected using a 128-element linear array

(M = 128). The transducer center frequency is 3.5 MHz and the sampling rate is

13.8889 MHz. The sound velocity is 1480 m/s. In Figure 4.8, we show the ultrasound

images obtained via DAS, SCB and RCB using only the central 32 and 64 transducers

of the array. This can be regarded as using smaller arrays with 32 and 64 elements,

respectively. For RCB, we choose e = 2. By using a 50 dB display dynamic range

in Figures 4.8(a)-4.8(f), we can discover the cysts from the images obtained via SCB

and RCB. However, the cysts are invisible from the images obtained via DAS (we see

instead some arcs). After checking carefully, we find out that, there are some samples

with very high amplitudes close to the end of several A-scan data sequences, which

results in very large DAS estimates at the arcs and make the cysts indistinguishable

from the background. As expected, the more transducers used, the more arcs will

appear in the DAS images. By increasing the display dynamic range to 80 dB for

DAS in Figures 4.8(g)-4.8(h), the cysts become visible. SCB and RCB do not have

such a problem since they can adjust their weight vectors and suppress the signals

corresponding to the abnormal samples in the A-scan data. Comparing the RCB

images and the SCB images, we notice that the cysts are much more discernible in

the former than in the latter. Hence RCB has better contrast than SCB.

We examine next a heart phantom. The data are collected using a 64-element

linear array (M = 64). The transducer center frequency is 3.333 MHz and the

sampling rate is 17.76 MHz. The sound velocity is 1540 m/s. In Figure 4.9, we show

the ultrasound images obtained via DAS, SCB and RCB using the entire array and









only the central 32 elements of the array, respectively. For RCB, we choose e = 2.

Indeed we cannot find the heart phantom at all from the poor SCB images. However,

the heart phantom is apparent in the DAS images and the RCB images since DAS

and RCB are much more robust against steering vector errors and small sample sizes

than SCB. It is clear that DAS has poor resolution and high sidelobes. Due to its

inherent adaptiveness, RCB does not suffer from these problems.

Finally, we consider in vivo imaging of a rat mammary tumor. The data are

collected using a 64-element linear array (M = 64). The transducer center frequency

is 2.6 MHz and the sampling rate is 25 MHz. The sound velocity is 1500 m/s.

In Figure 4.10, we show the ultrasound images obtained via DAS, SCB and RCB

using the entire array and only the central 32 elements of the array, respectively. As

expected, for each method applied to the same data, the more transducers used, the

better resolution we have. For DAS, we compare its images obtained using monostatic

data and multistatic data. Clearly, the images based on the multistatic data are

much better than those based from the monostatic data, as expected. Comparing

the images obtained using the same multistatic data, we see that the SCB images

are worse than those of DAS and RCB due to the fact that SCB is vulnerable to

the steering vector errors and small sample sizes. As can be seen, RCB has better

imaging results than DAS in terms of resolution and sidelobes. In particular, the

RCB image with M = 32 is very similar to the DAS image with M = 64, and it is

superior to the DAS image with M = 32.


































Depth Depth
(a) (b)


0 0


-4-4
6 5 10 15 20 5 10 15 20
Depth Depth
(c) (d)

Figure 4.1: Simulated ultrasound imaging results obtained via time-reversal based
methods. (a) MF, (b) SCB, (c) MUSIC, and (d) RCB with E = 1g112/100. A 30 dB
display dynamic range is used. There are no transducer position errors.






























2 2




5 10 15 20 5 10 15 20
Depth Depth
(a) (b)
0*0






-2 -2
-4 -4

5 10 15 20 5 10 15 20
Depth Depth

(c) (d)

Figure 4.2: Simulated ultrasound imaging results obtained via time-delay based meth-
ods. (a) DAS (monostatic data), (b) DAS (multistatic data), (c) SCB, and (d) RCB
with e = 1. A 30 dB display dynamic range is used. There are no transducer position
errors.






























S-2__-1 (__ +


5 10 15 20 5 10 15 20
Depth Depth
(a) (b)

424





-4-4

5 10 15 20 5 10 15 20
Depth Depth
(c) (d)

Figure 4.3: Simulated ultrasound imaging results obtained via time-reversal based
methods. (a) MF, (b) SCB, (c) MUSIC, and (d) RCB with e = IIg||2/10, A 30
dB display dynamic range is used. The transducer position errors are Gaussian
distributed with variance equal to 0.01A/2.





























Depth Depth
(a) (b)


Depth Depth
(c) (d)
Figure 4.4: Simulated ultrasound imaging results obtained via time-delay based meth-
ods. (a) DAS (monostatic data), (b) DAS (multistatic data), (c) SCB, and (d) RCB
with e = 1. A 30 dB display dynamic range is used. The transducer position errors
are Gaussian distributed with variance equal to 0.01A/2.



















60 E60







120 120
-50 0 50 -50 0 50
Lateral distance (mm) Lateral distance (mm)

(a) (b)

30 30
40 40
50 50
160
70 0
i 7o




110
120
-50 0 50 -40 -30 -20 -10 0 10 20 30 40
Lateral distance (mm) Lateral distance (mm)

(c) (d)

30 30
40 40
50 50



80 -


100 100
110 110
120 120
-50 0 50 -50 0 5o
Lateral distance (mm) Lateral distance (mm)
(e) (f)

Figure 4.5: Experimental ultrasound imaging results for the wire phantom. (a) DAS
with M = 64, (b) DAS with M = 128, (c) SCB with M = 64, (d) SCB with M = 128,
(e) RCB with M = 64, and (f) RCB with M = 128. A 50 dB display dynamic range
is used. For RCB, e = 10.








83







30 30
40 40
50 5
E E

0 d a n 70
a Wc
00


100 100
110 110
12 120
-50 0 50 -50 0 50
Lateral distance (mm) Lateral distance (mm)

(a) (b)


30 30
40 40
50 50
60 60
70 70



100 100
110 110
120 120
-50 0 50 -50 0 50
Lateral distance (mm) Lateral distance (mm)

(c) (d)

30 30
40 40












120 120
50 0 50 -50 0 50
















and c = 20. A 20 dB display dynamic range is used.







84















20 20

40 40

2C
S60 E 60








-60 -40 -20 0 20 40 60 -60 -40 -20 0 20 40 60
Lateral distance (mm) Lateral distance (mm)
(a) (b)






C C
80 o

100 100

120 120

-60 -40 -20 0 20 40 60 -60 -40 -20 0 20 40 60
Lateral distance (mm) Lateral distance (mm)

(c) (d)

Figure 4.7: Experimental ultrasound imaging results for the wire phantom with M =
64. (a) DAS (monostatic data), (b) DAS (multistatic data), (c) SCB, and (d) RCB
with c = 4. A 50 dB display dynamic range is used.
























50 50
E 60 60
70 70



100 00

120 20
-50 0 50 -50 0 50
Lateral distance (mm) Lateral distance (mm)
(a) (b)

30 30
40 40
50 50







110 10)

-50 0 50 -50 0 50
Lateral distance (mm) Lateral distance (mm)

(c) (d)

Figure 4.8: Experimental ultrasound imaging results for the cyst phantom. (a) DAS
with M = 32, (b) DAS with M = 64, (c) SCB with M = 32, (d) SCB with M = 64,
(e) RCB with M = 32, (f) RCB with M = 64, (g) DAS with M = 32, and (h) DAS
with M = 64. A 50 dB display dynamic range is used for (a)-(f), and a 80 dB display
dynamic range is used for (g)-(h). For RCB, e = 2.







































Lateral distance (mm)

(e)


u
Lateral distance (mm)

(f)


-50 0 50 -50 0
Lateral distance (mm) Lateral distance (mm)

(g) (h)


Figure 4.8: Continued.


50


a







87









6060
E E
80 80
1oo 100

140o 40
160 160

-80 -60 -40 -20 0 20 40 60 80 -80 -60 -40 -20 0 20 40 60 80
Lateral distance (mm) Lateral distance (mm)
(a) (b)

20 20

40 40






120 1 20

140 140



Lateral distance (mm) Lateral distance (mm)
(c) (d)

20 20

40 40




100 100

120 120
140 140

160 160
-80 -60 -40 -20 0 20 40 60 80 -80 -60 -40 -20 0 20 40 60 80
Lateral distance (mm) Lateral distance (mm)

(e) (f)

Figure 4.9: Experimental ultrasound imaging results for the heart phantom. (a) SCB
with M = 32, (b) SCB with M = 64, (c) DAS with M = 32, (d) DAS with M = 64,
(e) RCB with M = 32, and (f) RCB with M = 64. A 50 dB display dynamic range
is used. For RCB, e = 2.























30 30

E E
40 -40


50o C S
55 55
60 60
65 -
-20 -15 -10 -5 0 5 10 15 20 -20 -15 -10 -5 0 5 10 15 20
Lateral distance (mm) Lateral distance (mm)
(a) (b)

25 25
30 30


40 40
45 ._






-20 -15 -10 -5 0 5 10 15 20 -20 -15 -10 -5 0 5 10 15 20
Lateral distance (mm) Lateral distance (mm)
(c) (d)

Figure 4.10: Experimental ultrasound imaging results for the rat tumor. (a) DAS
(monostatic data) with M = 32, (b) DAS (monostatic data) with M = 64, (c) DAS
(multistatic data) with M = 32, (d) DAS (multistatic data) with M = 64, (e) SCB
with M = 32, (f) SCB with M = 64, (g) RCB with M = 32 and e = 2, and (h) RCB
with M = 64 and e = 4. A 50 dB display dynamic range is used.








89




















25

30 30

35 35
E E
40 40




55 55


65
-20 -15 -10 -5 0 5 10 15 20 -20 -15 -10 -5 0 5 10 15 20
Lateral distance (mm) Lateral distance (mm)

(e) (f)

25 25

30 30

35-35
E E




:N 50 50
55 55

60 60

65 65
-20 -15 -10 -5 0 5 10 15 20 -20 -15 -10 -5 0 5 10 15 20
Lateral distance (mm) Lateral distance (mm)

(g) (h)

Figure 4.10: Continued















CHAPTER 5
CONCLUSIONS AND FUTURE WORK
5.1 Conclusions

We summarize the issues addressed in the dissertation and the contributions

we made. We have developed a novel robust Capon beamformer (RCB) by com-

bining the covariance fitting formulation of the standard Capon beamformer (SCB)

with an ellipsoidal (including flat ellipsoidal) uncertainty set of the array steering

vector. We have shown how to efficiently compute RCB at a comparable compu-

tational cost with that associated with SCB. The data-adaptive RCB is much less

sensitive to steering vector mismatches than SCB and yet it can retain the appealing

properties of SCB including better resolution and much better interference rejection

capability than the standard data-independent beamformer. We have shown that the

RCB belongs to the class of diagonal loading approaches but the amount of diagonal

loading can be precisely calculated based on the uncertainty set of the steering vec-

tor. We have provided deep insights into the relationships among the recent three

robust adaptive beamformers. We have shown that, although the robust adaptive

beamformers obtained from the two different formulations appear to be rather differ-

ent, they are equivalent in terms of the signal-to-interference-plus-noise ratio (SINR).

Furthermore, we have demonstrated some unique features of RCB. Due to the fact

that both the power and the steering vector of the signal of interest are treated as

unknowns, there is a "scaling ambiguity" in the signal covariance term. This problem

is easily overcome in the RCB framework since the estimated steering vector can be

readily scaled. We have also compared these three beamformers for both the cases

of degenerate and non-degenerate ellipsoidal uncertainty sets of the steering vector,









and demonstrated the impressive advantages of RCB over the other two methods:

simpler implementation, lower computational complexity and more accurate power

estimation.

In addition to the theoretical study, we have applied the RCB algorithms to

various applications and extended RCB to accommodate the specific goals and re-

quirements. For the application of aero-acoustic imaging, we have devised a constant-

beamwidth RCB and a constant-powerwidth RCB for consistent acoustic imaging,

which means that for an acoustic wideband monopole source with a flat spectrum

the acoustic image for each frequency bin stays approximately the same. For the

application of ultrasound imaging, we have developed a time-delay based RCB and a

time-reversal based RCB, which have the desirable features including high resolution,

low sidelobes and robustness against steering vector errors. The effectiveness of RCB

and its various extensions has been demonstrated by simulated and experimental

examples.

5.2 Future Work

Several possible directions for future work are as follows.

The dissertation focuses on the batch mode (block adaptation) implementation

of RCB based on the stationary assumption. In the batch mode, a sample covariance

matrix needs to be constructed from a temporal block of array data before calculating

the RCB weight vector. In the nonstationary or time-varying situations such as the

target tracking applications, it is more suitable to implement RCB in the online mode

(continuous adaptation), which means that the weights are adjusted as the data are
sampled.

Throughout this dissertation, we have assumed a rank-one signal covariance

matrix. This assumption may be violated in sonar and wireless communications

where scattering sources or fluctuating wavefronts may exist. Hence it is of interest







92


to extend RCB to the general case of a signal subspace with a dimension larger than

one.

The performance of RCB will degrade in the presence of coherent signal and

interference. Preliminary work has been done to extend RCB for a simple case

where a priori knowledge of the directions of arrival of the pair of coherent signal

and interference is given. A more general case, where several coherent multipaths

impinge on the array from unknown directions of arrival with respect to that of the

signal, deserves attention for the future work.















APPENDIX A
LINKING RCB TO VOROBYOV AND COLLEAGUES' APPROACH

We repeat our optimization problem:


mina*R-la subject to Ia a|12 = e. (A.1)
a

Let ao denote the optimal solution of (A.1). Let

R-lao
wo la (A.2)
a*R-lao

We show below that the wo above is the optimal solution to the following SOCP

considered in [44]:

minw*Rw subject to w*a > /EIw|l + 1, Im (w*a) = 0. (A.3)
w

First we show that if I1111 <: /F, then there is no w that satisfies w*a >

V IlwII + 1. By using the Cauchy-Schwarz inequality, we have

V/EIIwII + 1 < w*a <_ VCIwI (A.4)

which is impossible. Hence the constraint in (2.29), which is needed for our RCB to

avoid the trivial solution, must also be satisfied by the approach in [44].

Next let

w = wo + y. (A.5)

We show below that the solution of (A.3) corresponds to y = 0.

Insertion of (A.5) in (A.3) gives:

2 1
min y*Ry + Re (y*ao) + I (A.6)
y aR-lao aR-lao




University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs