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 brothersinlaw 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 NonDegenerate 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 ConstantPowerwidth RCB . . . . . . . ..... ....... 44
3.4 ConstantBeamwidth 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 TimeDelay Based RCB ......................... 64
4.4 TimeReversal 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 dataindependent 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 dataindependent 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
signaltointerferenceplusnoise ratio (SINR). However, RCB also has some unique
features. By comparing these three beamformers for degenerate vs nondegenerate 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 aeroacoustic imaging, we devised a constant
beamwidth RCB and a constantpowerwidth 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 timedelay based RCB and a timereversal 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 informationbearing
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 spacetime 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 farfield data model to show some concepts of array signal
processing. Wideband signals in the nearfield 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 farfield.
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 farfield. We assume that the data and weights are complex
valued, since in many applications a quadrature receiver is used to generate inphase
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 Jdimensional
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)eiw ... HM(wc)eicrM]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 eiwcr2 .. eiwTM]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 eiW ... ei(M1)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 halfwavelength 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 delayandsum 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 delayandsum 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 dataindependent beamformer selects its weight vector
regardless of the incoming data. On the other hand, an adaptive (datadependent)
beamformer takes advantage of the incoming data and adjusts its weight vector as
a function of the data. Standard dataindependent beamformers include the delay
andsum approach and methods based on various dataindependent 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 datadependent Capon beamformer adaptively selects the weight vector
to minimize the array output power, subject to the linear constraint that the signalof
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 signaltointerferenceplusnoise ratio (SINR).
Then the performance of the Capon beamformer may become worse than that of the
standard dataindependent 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 smallsample 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: nondegenerate
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 aeroacoustic
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
nondegenerate 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 socalled 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 signalofinterest (SOI), and that the remaining
rankone 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 interferenceplus
noise vector for the nth snapshot. In the radar application, the SOIfree data vectors
are available, so the interferenceplusnoise 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 interferenceplusnoise 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 nondegenerate 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 timeconsuming
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:
Rlao
wo = Rl (2.16)
a*Rlao
Using (2.16) in Step (b) yields the following estimate:
a2 = a0 (2.17)
0 aORlao
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
semidefinite. The previous claim follows from the following readily verified equiva
lences (here R1/2 is the Hermitian square root of R'):
R a02aa* > 0
I a2R1/2aoaR12 > 0 ::
1 a2aRlao > 0 <>
o2 < aR a = 2 (2.19)
Sa;Rla
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
semidefinite.
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 NonDegenerate Ellipsoidal Uncertainty Set
When the uncertainty set of the steering vector a is a nondegenerate 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)* C1 (a a) < 1 (2.20)
where a and C are given.
The RCB problem in (2.20) can be readily reformulated as a semidefinite
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 socalled 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*Rla subject to (a a)* C1 (a a) < 1 (2.22)
a
To exclude the trivial solution a = 0 to (2.20), we assume that
a*Ca > 1 (2.23)
Note that we can decompose any matrix C > 0 in the form:
1
C1= 1D*D (2.24)
E
where for some e > 0,
D = vC1/2 (2.25)
Let
S= Da, =Da, R =DRD* (2.26)
Then (2.22) becomes
mina*lRi 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*Rla subject to Ila a112 < e (2.28)
a
To exclude the trivial solution a = 0 to (2.28), we now need to assume that
lall2 > 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*Rla + v (la a12 e) (2.30)
where v > 0 is the realvalued 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*Rla 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*R1 + 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
(R1 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*Ra for any a E S (2.35)
Maximization of h2(v) with respect to v gives
(V +l (2.36)
which indeed satisfies
Ilo a12 = (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 eigendecomposition 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 lao 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 eigendecomposition 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 (v2I + 2vFl + r2)1 U*a
where the inverse of v2I + 2vlr + 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:
Rlfo
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 realvalued 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 signaltointerferenceplusnoise 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
*Cli^ > 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 soobtained 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 eigendecomposition 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 secondorder
cone program (SOCP) approach in [44]. The approach in [45] can be implemented
recursively by updating the eigendecomposition 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)*R1(Bu + a) subject to Ilull < 1 (2.53)
U
Note that
(Bu + a)*Rl(Bu + a) = u*B*R1Bu + a*R1Bu + u*B*R1a + a*R'a (2.54)
Let
R = B*R1B > 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 MoorePenrose pseudoinverse 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 IlR1all < 1, then the unique solution in (2.63) with 0 = 0, which is it = Rl',
solves (2.57). If Ili1'AI > 1, then f > 0 is determined by solving
2( I (ft PI)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) = IRil1I > 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 eigendecomposition of R (see (2.65)).
Step 3: If IIR1ll < 1, then set 0 = 0. If 11R111 > 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
R1 and the eigendecomposition of R. If L << M, then the complexity is mainly
due to computing R1. 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 signaltointerferenceplusnoise 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 nondegenerate 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
0adRldo
(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 = R1B(R+ PI) + Rla
= RB(B*R1B + PI) B*R1a + Ra
= (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 halfwavelength 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 MonteCarlo simulations are given. However, the beampatterns shown
are obtained using R from one MonteCarlo 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
Ia(0o) al2 < 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 oo 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/nsojl 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
S80 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
S40 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 ooo 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 oo 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
10 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 delayandsum
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 zeromean 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 dataindependent beamformer using
the assumed array steering vector as the weight vector. This approach is referred to
as the delayandsum 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 Delayandsum  Delayandsum
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 delayand
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 ra &~ 
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 +A3)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 . oo
a 15 15
S10 10
,,, l2 5 .. .''''"
5 /
10 A RCB (Flat ellipsoid) 10o RCB (Flat ellipsoid)
 RCB (Sphere) w RCB (Sphere)
So 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 wellknown that as
the frequency increases, the beamwidth of both dataindependent and dataadaptive
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 constantbeamwidth 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
constantbeamwidth 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 delayandsum (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 constantbeamwidth
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 frequencydependent uncertainty parameter for the array steering vector;
we refer to the soobtained beamformer as the constantpowerwidth 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 Ipoint
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 complexvalued 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 farfield 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 rootmeansquared 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 timedomain sequence
with I samples, and let Y(i), i = 1, 2, .. I, denote the frequencydomain sequence
obtained by applying an Ipoint FFT to {y(1)}. Let yi(l),l = 1,2,...,I, be the
narrowband timedomain sequence corresponding to the ith frequency bin of the
above frequencydomain 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 rootmeansquared 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 rootmeansquared 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 constantbeamwidth and
constantpowerwidth 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 interferenceplusnoise covari
ance matrix. If xg = xo and the signaltointerferenceplusnoiseratio (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 twodimensional (2D) array
imaging, in which the beamwidth or powerwidt is defined as the diameter of a circle
having the same area as the 3dB contour of the main lobe of a 2D 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 ConstantPowerwidth 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 signaltonoiseratio (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 frequencydependent 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 ConstantBeamwidth 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 constantbeamwidth 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 semidefinite 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 dataindependent 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 
00 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 wavenumberlength product kD,, where k is the wavenumber and D, is the
diagonal distance between the elements of the nth cluster. The wavenumberlength
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
xy plane, with center location at (0, 0, 0). Note that inch is used as the unit for the
3D coordinates.
We assume that the distance between the array and the source is known and
plot the 2D 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 8192point FFT on the nonoverlapping 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 3dB 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 3dB 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 
88
?12\ ?12 .
"). . 0)
..... . ._ ._ . .
4 4L
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 3dB 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 3dB 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 130
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 dataindependent 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 dataindependent delayandsum (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 dataindependent
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 pulseecho 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
signaltonoise ratios of the individual narrowband beamformers [92], [93]. To make
better use of the large bandwidth, several methods have been proposed for wideband
farfield processing using passive arrays. For example, the coherent signalsubspace
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 spacetime 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 nearfield 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
pseudocovariance 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 reemitted 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 pseudocovariance matrix constructed in the spacetime
domain similarly to [96], while the timereversal 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 timedelay based RCB can tolerate the misalignment
of data samples and the timereversal 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 pulseecho mode to
probe the unknown propagating medium, i.e., pulses are emitted and data sequences
of the backscattered 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 Ascan data for all possible
transmitter and receiver pairs of the array. Specifically, the Ascan data sequence
Pi,(t) is obtained by recording the backscattered echo at the jth transducer, after
transmitting a pulse from the ith transducer. Assume that each Ascan 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 Ascan 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 onaxis 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 (offaxis 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 pseudocovariance 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 spacetime matrix R(p) as a pseudocovariance
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 timedelay based SCB is to pass the onaxis signal
(the signal on the focal point) without distortion, and at the same time minimize the
contributions from the offaxis signals (other signals on the arcs). However, in practice
the timedelay 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 smallsample 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 onaxis or offaxis. It cannot adapt to the incoming
data and hence has poor resolution and offaxis signal suppression capability.
From the above discussions, we see that the timedelay 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 (kiy 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 timereversal based adaptive beamformer may degrade severely. This motivates us
to consider a robust adaptive beamformer. Similar to the timedelay based method,
we can apply RCB by allowing the Green's vector to be within an uncertainty set as
in (4.8).
4.3 TimeDelay Based RCB
In this section we focus on the timedelay based RCB, with additional dis
cussions on timedelay 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 timedelay 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:
R1a R1Mxl
S= = (4.14)
a*Rla l1xlR11Mxl
Using (4.14) in Step (b) above yields the following estimate of a2:
2 = 1 1(4.15)
a*Rla 1MxlR11Mx1
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 signalofinterest (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*Rla subject to Ila a2 < (4.17)
a
To exclude the trivial solution a = 0 to (4.17), we assume that
IIa 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 reformulate (4.17) as the following quadratic problem
with a quadratic equality constraint:
mina*Rla 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 lefthand 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):
RlA
w = (4.24)
a*R1f
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 eigendecomposition). Therefore, the computational
complexity of RCB is comparable to that of SCB.
When the pseudocovariance matrix R in (4.13) is replaced with an identity
matrix corresponding to the pseudocovariance matrix for spatially white noise, it
becomes a delayandsum 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 dataindependent 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 TimeReversal Based RCB
In this section we discuss the timereversal based RCB and several other al
gorithms including the timereversal based SCB, matched field (MF), and multiple
signal classification (MUSIC).
We consider herein a coherent approach, which exploits the timebandwidth
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 timereversal based methods, whose formula
tions are analogous to the timedelay based counterparts in the previous section. To
simplify the notation, we omit the argument wo in T(wo) and g(p, wo).
The timereversal based RCB has the form:
mina*Tla 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 IalI = llg(p)ll as follows:
a = II(p)lla/lal (4.38)
The remaining steps of the timereversal based RCB are in direct analogy to those
of the timedelay based RCB, and are hence omitted herein for brevity.
The timereversal 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)TIg(p)
The timereversal 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 timereversal 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 timedelay based
DAS in (4.25).
Once we have obtained the weight vector w for the timereversal based SCB,
RCB and MF, we can estimate the magnitude square of the onaxis 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 timereversal based imaging. MUSIC
was a superresolution 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 wellknown, 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 timereversal and timedelay 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 intersensor 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 freespace 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 timedelay based
methods, each Ascan data sequence is preprocessed by a matched filter based on
the probing pulse. For the timereversal 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 timereversal and
timedelay based methods, respectively, for the case of no transducer position errors.
We see from Figure 4.1 that the timereversal based MF has poor resolution and high
sidelobes. Although the timereversal based SCB image has peaks corresponding to
the targets, it also contains some sidelobes like that of MF. Both MUSIC and the
timereversal based RCB have good crossrange 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 timereversal based RCB outperforms the other
three methods. From Figure 4.2, we see that the timedelay 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
timedelay based RCB has the best imaging result among all the timedelay based
methods considered here. As can be seen from Figures 4.1 and 4.2, the timedelay
based approaches have better range resolution but worse crossrange resolution than
the timereversal 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 zeromean 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 timedelay based methods
and uncertainty in the Green's function for the timereversal based methods. Hence
there are array steering vector errors. As a result, degradations occur in the ultra
sound images. Note that the timereversal based RCB and the timedelay 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 timedelay 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 UrbanaChampaign
(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 128element 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 64element 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 Ascans 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 128element 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 Ascan 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 Ascan 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 64element
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 64element 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
44
6 5 10 15 20 5 10 15 20
Depth Depth
(c) (d)
Figure 4.1: Simulated ultrasound imaging results obtained via timereversal 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 timedelay 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.
S2__1 (__ +
5 10 15 20 5 10 15 20
Depth Depth
(a) (b)
424
44
5 10 15 20 5 10 15 20
Depth Depth
(c) (d)
Figure 4.3: Simulated ultrasound imaging results obtained via timereversal based
methods. (a) MF, (b) SCB, (c) MUSIC, and (d) RCB with e = IIg2/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 timedelay 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
3535
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 dataadaptive 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 dataindependent 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 signaltointerferenceplusnoise 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 nondegenerate 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 aeroacoustic imaging, we have devised a constant
beamwidth RCB and a constantpowerwidth 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 timedelay based RCB and a
timereversal 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 timevarying 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 rankone 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*Rla subject to Ia a12 = e. (A.1)
a
Let ao denote the optimal solution of (A.1). Let
Rlao
wo la (A.2)
a*Rlao
We show below that the wo above is the optimal solution to the following SOCP
considered in [44]:
minw*Rw subject to w*a > /EIwl + 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 CauchySchwarz 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 aRlao aRlao
