NAVAL POSTGRADUATE SCHOOL
Monterey, California
THESIS
DOA ESTIMATION BY EIGENDECOM POSITION
USING SINGLE VECTOR
LANCZOS ALGORITHM
by
Daniel E. Gear
June 1989
Thesis Advisor Murali Tummala
Approved for public release; distribution is unlimited.
T244080
classified
; Report Security Classification Unclassified
REPORT DOCUMENT ATION PAGE
Restrictive Markings
Security Classificaiio
Declassification Downgrading Schedule
Distribution Availability of Report
Approved for public release; distribution is unlimited.
'erforming Organization Report
1 Report Number(s)
Name of Performing Orgamzati
aval Postgraduate School
6b Office Symbol
i if apfluable) 62
7a Name of Monitoring Orgamzati
Naval Postgraduate School
Address icily, state, and ZIP code)
[onterey. CA 939435000
7b Address (city, state, and ZIP code)
Monterey, CA 939435000
Name of Funding Sponsoring Or
Identification Number
Address (city, state, and ZIP code)
e of Funding Numbers
Program Element No  Project No [Task No  Work Unit Accession No
Title ( include security classificaiio
ANCZOS ALGORITHM
DOA ESTIMATION BY EIGENDECOMPOSITION USING SINGLE VECTOR
k Personal Authon's) Daniel E. Gear
la Type of Report
paster's Thesis
14 Date of Report (
June 1989
) Supplementary Notation The views expressed in this thesis are those of the author and do not reflect the official policy or po
tion of the Department of Defense or the U.S. Government.
IS Subject Terms (continue on reverse if nereis an and identify by Mod
direction of arrival, Lanczos, eigendecomposition.
1 Abstract i continue on reverse if necessary and identify by block number)
Subspace methods of solving spectral estimation and direction of arrival (DOA) problems involve finding the eigenvalues
nd eigenvectors of correlation matrices. Using the Lanczos algorithm some of the extreme eigenvalues and eigenvectors can
>e approximated without requiring the entire matrix decomposition theoretically saving many computations.
This thesis develops a model and a form of the Lanczos algorithm to solve the DOA problem. The relationship of the
lumber of eigenvectors used to the accuracy of the results in a low signal to noise ratio example are examined.
!0 Distribution Availability of Abstract
X unclassified unlimited Q same a:
□ DTIC use
21 Abstract Security Classification
Unclassified
>2a Name of Responsible Individual
vlurali Tummalu
22b Telephone i include Art
(408) 6462645
22c Office Symbol
62Tu
>D FORM 1473,84 MAR
i of this page
Unclassified
Approved for public release; distribution is unlimited.
DOA ESTIMATION BY EI GEN DECOMPOSITION
USING SINGLE VECTOR
LANCZOS ALGORITHM
by
Daniel E. pear
Lieutenant Commander, United States Navy
A.B., University of North Carolina, 1977
Submitted in partial fulfillment of the
requirements for the degree of
MASTER OF SCIENCE IN ELECTRICAL ENGINEERING
from the
NAVAL POSTGRADUATE SCHOOL
JUNE 1989
ABSTRACT
Subspace methods of solving spectral estimation and direction of arrival (DOA)
problems involve finding the eigenvalues and eigenvectors of correlation matrices. Using
the Lanczos algorithm some of the extreme eigenvalues and eigenvectors can be ap
proximated without requiring the entire matrix decomposition theoretically saving many
computations.
This thesis develops a model and a form of the Lanczos algorithm to solve the DOA
problem. The relationship of the number of eigenvectors used to the accuracy of the
results in a low signal to noise ratio example are examined.
C.l
TABLE OF CONTENTS
I. INTRODUCTION 1
II. DIRECTION OF ARRIVAL 3
A. SPECTRAL ESTIMATION 3
B. BEAMFORMING 8
C. SUBSPACE METHODS 13
1. Introduction 13
2. MUSIC 15
3. ESPRIT 18
4. Other Subspace Methods 21
III. THE LANCZOS ALGORITHM 22
A. CLASSICAL LANCZOS AND ITS VARIANTS 23
B. IMPLEMENTATION 27
1. Single Vector Lanczos Algorithm 28
2. Other Methods 39
IV. RESULTS 42
A. APPROACI I 42
B. EXPERIMENT SET UP 43
V. CONCLUSIONS AND RECOMMENDATIONS 67
LIST OF REFERENCES 68
INITIAL DISTRIBUTION LIST 73
LIST OF FIGURES
Figure 1. Wavefront 9
Figure 2. Simple delayandsum beamformer 10
Figure 3. Spatial Frequency 12
Figure 4. Steering vector for 3 sensors and 1 signal 16
Figure 5. Signal Subspace for 2 signals 17
Figure 6. PSD of first eigenvector 32
Figure 7. PSD of second eigenvector 33
Figure 8. Overlayed PSDs of first 5 eigenvector 34
Figure 9. Eigenvector averaging 35
Figure 10. Spectral averaging 35
Figure 1 1. Spectral product for 2 eigenvectors 36
Figure 12. Spectral product for 5 eigenvectors 36
Figure 13. Spectral product for 5 eigenvectors, dB 37
Figure 14. Spectral product for 5 eigenvectors, 5 dB 37
Figure 15. Spectral product for 5 eigenvectors, 5 dB 38
Figure 16. Spectral product for 10 eigenvectors, 5 dB 38
Figure 17. Spectral product for 15 eigenvectors, 5 dB 39
Figure IS. A physical implementation 44
Figure 19. Case 1 5 dB, 1 target at 18 ° 45
Figure 20. Case 1 5 dB. 1 target at 81 ° 46
Figure 21. Case 1 5 dB, 2 targets at 36 ° and 3S ° 47
Figure 22. Case 1 5 dB, 3 targets 4S
Figure 23. Case 2 dB, 1 target at 18 ° 50
Figure 24. Case 2 dB, 1 target at 81 ° 51
Figure 25. Case 2 dB, 1 target at 81 ° 52
Figure 26. Case 2 dB, 2 targets at 36 ° and 38 ° 53
Figure 27. Case 2 dB, 3 targets at °, 36 ° and 88.2 ° 54
Figure 28. Case 3 5 dB, 1 target at 18 ° 55
Figure 29. Case 3 5 dB, 1 target at 81 ° 56
Figure 30. Case 3 5 dB, 2 targets at 36 ° and 38 ° 57
Fisure 31. Case 3 5 dB, 2 tarsets at 36 ° and 40 ° 58
Figure 32. Case 3 5 dB, 3 targets at °, 36 ° and 88.2 ° 59
Figure 33. Case 4 5 dB, 3 targets at IS °, 36 ° and 41.4 ° 60
Figure 34. Case 4 5 dB, 3 targets at 18 °, 36 ° and 41.4 ° 61
Figure 35. Case 4 5 dB, 3 targets at 18 °, 36 ° and 41.4 ° 62
Figure 36. Case 4 5 dB, 3 targets at 18 °, 36 ° and 41.4 ° 63
Figure 37. Case 4 5 dB, 3 targets at 18 °, 36 ° and 41.4 ° 64
Figure 38. Case 4 5 dB, 3 targets at 18 °, 36 ° and 41.4 ° 65
Figure 39. Case 4 5 dB, 3 targets at 18 °, 36 ° and 41.4 ° 66
I. INTRODUCTION
Recently, a number of signal processing techniques have been developed that in
volve resolving a received signal into socalled signal and noise subspaces. The ability
to perform spectral estimation or direction of arrival from the determination of the noise
subspace has been shown in the works of Pisarenko, Schmidt, Kay, Kailath, and others.
The methods van' in the manner in which the subspace is reached and how the resulting
eigenvectors are calculated. [References 1,2, and 3].
To achieve better results (detection at lower signaltonoise ratios, better resolution,
finer accuracy, less bias), more samples are required. This leads to larger, more complex
matrices that better describe the signal and noise subspaces. Determination of the sub
spaces requires an eigendecomposition of an estimated correlation matrix. The general
eigendecomposition of a matrix requires computations 0(n 2 ), thus the larger matrices
require many more operations. Once the matrix is decomposed into its eigenvalues and
eigenvectors the proper set of eigenvectors must be used to find the resulting spectrum.
Hence, estimation is required to split the eigenvalues into signal and noise subsets.
The difficulties in the procedures can be stated with two questions.
• Where is the threshold between signal and noise eigenvalues (and thus which, or
how many eigenvectors are used)?
• What method should be used to find those eigenvectors in a timely fashion?
Proposed here is a method which will accurately estimate some of the eigenvalues
and eigenvectors of a matrix without performing the entire decomposition. Computa
tional savings are realized when only a partial decomposition is required. The
eigenvectors used correspond to the minimum eigenvalues of the matrix. With an
overspecified matrix (dimension much larger than the number of signals present), these
minimum eigenvalues will all correspond to the noise subspace. Each noise eigenvector
contains all the information to find the actual spectrum, although spurious results will
also be apparent (because the matrix is over specified). Several methods of handling
these spurious peaks are illustrated, including eigenvector averaging, spectral averaging
and using the spectral product.
This thesis is organized in five chapters. Chapter II, Direction of Arrival, discusses
classical spectral estimation theory and how it applies to the linear array problem
(beamforming). Subspace methods starting with Pisarenko Harmonic Decomposition
and proceeding to MUSIC and ESPRIT are discussed in detail. Chapter III, The
Lanczos Algorithm, includes all the basic theory required to describe this
eigendecomposition method. There we also compare several methods to negate the
spurious peaks with the proposed spectral product giving the best results. Results of the
algorithm for numerous cases are given in Chapter IV. The last chapter summarizes the
results and advantages of this Lanczos subspace method. This chapter also includes
some recommendations for future work and lists some possible applications.
II. DIRECTION OF ARRIVAL
Direction of arrival is a form of spectral analysis, performing frequency detection
and resolution in the spatial domain vice the conventionally considered time domain.
The signals incident on an array are analyzed, and, if the presumptions of the analysis
are valid, the correct bearings to the sources are determined. Formerly it was only pos
sible to analyze the output of an array by conventional Fourier techniques. More re
cently numerous methods have been developed which enhance the ability to accurately
determine the spectral and angular resolution.
This chapter summarizes some of the salient points of spectral estimation before
discussing classical direction of arrival array processing (Bartlett beamforming).
Projectiontype superresolution subspace methods are then discussed, starting with
Pisarenkos Harmonic Decomposition. MUSIC and ESPRIT are discussed in detail and
several other subspace methods are mentioned.
A. SPECTRAL ESTIMATION
Spectral estimation is the term used to describe efforts to find the frequency com
ponents of a signal sampled in time. Two conditions that are required for the remainder
of this thesis are that the processes considered are wide sense stationary (WSS) and
ergodic. The assumption of WSS means that the process has finite mean and that the
autocorrelation is a function of the distance, or lag, between two samples and not of the
position of the samples themselves. Ergodicity allows time averages to be used to de
termine ensemble averages.
The autocorrelation function of the signal x(t) is
R xx (m) = <x(n + m),x(n)> = ^l 2N \ { / x(n + m)x{n)l (1)
where x(n) are the individual samples of the signal. When only a finite number of sam
ples, M, are taken, the above definition must be modified. The sample autocorrelation
function is defined as
««(*)
Mlk
+ k)x{n), for 0<*<(A/1)
(2)
It is easily shown that [Ref. 4: pp. 5658]
R^k) = ^ A), for  (M  1) < k <
and
*«(<>) > *«(*). for all k
The sample autocorrelation matrix R„ is
R^O) R xx {\)
*«(1) 4x(0)
K^O K^Af1)
R xx (M+\)
««(0)
(3)
(4)
The power spectral density (PSD) is a measure of power per unit frequency. A plot
of PSD versus frequency will show the relative power at all the frequencies present. It
also describes the properties of the noise in the signal, i.e., white noise will have a flat
PSD showing that all frequencies are equally represented [Ref. 1: pp. 53]. The PSD is
given by
SJJ)
+2
R xx (m)e 
where T is the sampling period.
The periodogram method of estimating the true PSD is one of the earliest and most
widely used [Ref. 5: pp. 58]. The periodogram is defined as the ztransform of the
sample autocorrelation matrix evaluated on the unit circle [Ref. 4: pp. 5253]
*=(Afl)
A/1
z
(6)
It may also be found by directly ztransforming the original data sequence x(n)
S xx (z) = jj X(z)X(z ] ) where X(z)
The periodogram spectrum is found by substituting z = e i2nfr ,
M
(7)
Lif) = jr\X(f)\ 2 = 7
I
x(n)e •
(8)
Data is often run through a computationally efficient Fast Fourier Transform (FFT) to
find the periodogram spectrum.
The measures of effectiveness of a spectral estimator are based on comparisons of
"resolution, detectability, bias and variance. Resolution relates to the ability of the esti
mator to separate two separate spectral lines that are closely spaced in frequency. The
capacity to locate a low energy signal is measured in detectibility. Bias is the tendency
of the observed peak to be shifted off the true spectral line. Variance is a measure of the
propensity to keep the true spectral shape, or how quickly the PSD falls of! away from
the true peak. Different spectral estimators are sometimes compared in terms of
robustness, or ability to function well with various types of signals and noise.
If individual signals are narrowband the resolution criterion is said to be achieved
when direct observation of the spectrum leads to the correct determination of the num
ber of signals. Separate peaks are not required to determine that two signals are present
[Ref. 6]. Resolution is inversely proportional to the length of the data samples. With/
as the sampling rate and T the sample period, or time between samples, T = ~, the fre
quency resolution is given by
/,'
V  17  iff (9 »
Thus, with the periodogram, better resolution can only be achieved with longer record
lengths.
The above description is based on no windowing (or rectangular windowing) of the
data samples. The windowing in the time domain is seen as convolution in the frequency
domain. The rectangular windows, for example, transform into sine functions in the
frequency domain resulting in high sidelobe distortion known as leakage in the frequency
domain. High sidelobes result in many false peaks, making it difficult to discern the true
spectral peaks, so detectibility suffers.
Nonrectangular windows are used to taper the data to minimize the discontinuities
that cause the high sidelobes. Common windows include the Bartlett and Hamming.
The lower sidelobes come at a cost of resolution, so two or more signals close to each
other in frequency may merge into one. Resolution may be boosted by lengthening the
sample time, but the increased record length may violate the requirement for
stationarity. More can be found on the subject of windows in References 7 and 8.
Because of the difficulty in meeting the requirements for both detectibility and re
solution, parametric methods have been developed that produce increased resolution
with short data lengths. The parametric methods are based on determining an appro
priate model for the process that produced the data. If the process can be effectively
modeled, then more reasonable assumptions may be made about the data's behavior
outside of the window (i.e., these data points do not have to be set to zero). Once the
appropriate model is chosen to represent the process, the parameters are estimated from
the available data, inserted into the model, and the power spectral density expression for
the respective model is determined. A few common parametric methods are summarized
below. [Ref 1: pp. 10610S, Rcf5: pp. 172174]
Many random processes that are encountered in practice may be represented by the
linear difference equation
p <?
• \ a(k)x(n  k) + y
x(n) =  > a(k)x(nk) + > b(k)u{n  k) (10)
where u{n) is the input driving sequence and x{ri) is the output sequence. The a(k) pa
rameters are the autoregressive coefficients and the b(k) are the moving average coeffi
cients. Equation 10 is thus known as the autoregressivemoving average (ARMA) model
or AKM.A r p,q) process and is the most general model with a rational transfer function.
The power spectral density that results from the ARMA model exhibits both sharp
peaks and deep nulls. If the autoregressive parameters of equation 10 are all set to zero
except a(0) = 1, the resultant model is the moving average (MA) process of order q. The
MA spectrum will have deep nulls, but relatively broad peaks. With the b(k) coeffi
cients of equation 10 set to zero except b(0) = /, the autoregressive (AR) process of order
p results. The AR spectrum will exhibit sharp peaks but will have broad nulls. [Ref 5:
pp. 1741 SI] Each of the models will exhibit greater accuracy and flexibility as the order
is increased. With a high enough order the AR method can approximate an ARMA or
MA process, and, likewise, a very high order MA model can be used to approximate an
ARMA or AR process. But if the model order is too high for the actual process, esti
mation errors will lead to nonzero coefficients and spurious peaks will result. Thus
proper estimation of model order is important [Ref 1: pp. 112113, pp. 19S207].
The parameters o[ these three models may be estimated by using the YuleWalker
equations utilizing the autocorrelation matrix of the available data stream [Ref 1: pp.
115118]
R XA a = r (11)
While the true autocorrelation function is impossible to determine, the maximum
likelihood estimator (ML) is one means to find good approximations of the parameters
for the AR model. The ML method uses a suitable estimate of the autocorrelation or
covariance matrix and then solves [Ref 1: pp. 1S5190J
Ca = c (12)
for the parameters a where C is the symmetric covariance matrix and c is an
autocorrelated vector.
The Burg method (maximum entropy) estimates reflection coefficients from the
process and then uses the Levinson recursion to find the estimated parameters [Ref 1:
pp. 228231].
Generally, the various AR spectral estimators except the autocorrelation method are
unbiased and have similar variance. The covariance and Burg methods are restricted to
ranges that keep the summation within the available data and do not assume zero pads
outside of the samples, thus taking advantage of the basis which led to the creation of
the parametric methods in the first place. [Ref 1: pp. 240252]
B. BEAMFORMING
A conventional approach to the direction of arrival (DOA) problem is through the
classical beamforming (Bartlett) technique. Basically, this is a measure of coherency of
the signals arriving at an array of sensors. The characteristics of an array are described
in terms of array gain, directivity, resolution, beamwidth, and sidelobes. These are based
on array size (number of sensors), sensor sensitivity, intersensor spacing, and post re
ception processing.
Assuming a narrowband point source in the far field, a plane wave will pass through
a linear array in a specified order, where the magnitude of the excitation on any indi
vidual sensor will be related to the phase of the signal at the instant of sampling. The
individual sensors will see different instantaneous magnitudes based on this phase dif
ference which is a function of the frequency (or wavelength), the DOA, and the spacing
between sensors. The difference in phase for two successive sensors for a linear array
with equally spaced sensors is
A0 = ~^sin0 (13)
where d is the distance between sensors. / is the signal wavelength, and 6 is the angle
between the wavefront and the array axis. This phase difference A</> is also known as the
normalized wavenumber, k [Ref. 4: pp. 341343], Figure 1 illustrates the array
wavefront interaction.
The output of a simple narrowband delayandsum beamformer is
M
I
e{t) = > x m {tx m ) (14)
where x„(t), m— 1, 2, ... , M is the signal at the mth sensor. The time lag, r mt between
two adjacent sensors is to make up for the propagation delay caused by the wavefront
having to travel the extra distance ds'mO. One can adjust the magnitude of the output
to plane waves arriving from a particular direction by introducing appropriate time
delays at each sensor prior to summing. This is known as "steering the array" or
"steering the beam". An illustration of a typical beamformer arrangement is shown in
Figure 2.
WAVEFRONT
ARRAY
Figure 1. Wavefront
In the multiple source case, especially if the undesired signals are interfering with the
detection of other sources, this technique may be modified to steer nulls instead of
beams thus minimizing output from "jammers". Nulls may be directed toward a single
known source to minimize the strength of the signal with respect to that source, with the
output being analyzed to determine what other sources are present. Another imple
mentation is to steer the nulls to minimize the total output. The analysis of the delays
will indicate the directions of multiple targets. [Rcf. 4: pp. 351355, Ref. 9 ]
Beamforming may also be analyzed in the frequency domain by transforming
equation 14 into the frequency domain
*(/)
I
*M*~
(15)
O O Sensors O
[ T ll H S Delays Q
\U!
OUTPUT
Figure 2. Simple delayandsum beamfoi nier
In vector notation we can write e = \\ T x where w are the weights and x are the outputs
of M sensors
where w =
and x =
* 2 W
,***
*m(/)
(16)
The steering vector s k is the phase relationship between the angle and the normalized
wavenumber k given in equation 13
and k
2nd
sin 9
(17)
It can be shown that the steering vector from an array with weights w, directed toward
an arbitrary direction 6, can be expressed in terms of the steering vector as a = s A [Ref.
4: pp. 343344].
Frequency domain beamforming is directly analogous to spectral estimation
decribed above. Spatial sampling has the requirement that sensor distances d must be
less than ).\2 apart to prevent "grating lobes" (or aliasing) corresponding to the Nyquist
rate in the frequency domain of f m3X <fJ2 . Longer arrays, containing more sensors,
will give better resolution and smoothing. This is similar to frequency resolution being
proportional to the data record length (A/a ~Trf)> t Re ^ 4: PP 341345. Ref. 10: pp.
48, Ref. 11 : pp. 293299]
The DOA is actually a relative comparison of observed spatial frequency and known
signal frequency. The spatial frequency is greatest on endfire, when the wavefront is
perpendicular to the array (6 = 90° or — ). Here the phase difference between adjacent
sensors is at a maximum. A simultaneous sampling of all sensors at one instant in time,
or snapshot, will indicate the maximum spatial frequency. Observation of the spatial
wave over time (with a known temporal sample rate) will indicate the end of the array
where the source is located.
On broadside (6 = or n), the wavefront excites each sensor identically. Spatial
sampling indicates no phase difference along the entire array, except for the effects of
additive noise, resulting in the computation of zero spatial frequency. Unfortunately,
linear arrays have an inherent ambiguity with broadside signals. The side of the array
at which the source is located cannot be determined without further information. Spatial
frequency is illustrated in Figure 3.
An extra requirement in the standard DOA problem is a priori knowledge of in
coming signal frequencies. Typically, this is handled via a bank of bandpass filters on
the output of the sensors. Data streams from the sensors at each desired center fre
ARRAY SENSOR
ARRAY SENSOR
Figure 3. Spatial Frequency: Two signals with SNR = 2dB,
0, = 1° and 2 = A5° for two snapshots at time = (a) /„ , (7>) /,. Note
the variation in 'DC level' as the snapshots sample the nearly broadside
signal at different phases.
quency (/J, f 2 , ...) are processed in parallel, resulting in spatial frequencies for each
temporal frequency. The DOAs are calculated by comparing these spatial frequencies
with the center frequencies of their respective filter banks.
Improvement in beamforming may be seen through the use of windows, weighing
each sensor output by the appropriate amount to narrow the beamwidth or lower
sidelobes, but, as discussed earlier, at a cost of lowering overall resolution.
C. SUBSPACE METHODS
1. Introduction
The projectiontype subspace method utilizes the properties of the
autocorrelation, covariance, or modified covariance data matrices and their
eigenvalue/eigenvector decomposition into signal components and noise components in
estimating the DOA. Generally, subspace methods use an assumed property of the data
to provide good resolution if the data fits the particular assumptions, even for extremely
short data sets. If the data (and signals) do not meet the assumptions, the results can
be quite misleading. The assumptions here call for white noise and a signal whose esti
mated autocorrelation matrix converges to the true autocorrelation matrix as the order
is increased.
For p complex sinusoids in additive complex white noise the combined
autocorrelation function of the signal plus noise is given by
R xx { fn ) = \ p^ff" + a\d{m) (18)
2*
We define R„ as the signal autocorrelation matrix of rank p, R Z2 = a\\ of full rank M
and give the autocorrelation matrix as
p
R xx = \ Pfitf + a 2 v l = R ss + R 22
(19)
where I is an A/by M identity matrix, R xx is of full rank M due to the noise contribution,
P, is the power of the /th complex sinusoid, a] is the noise variance, and
e, = [ 1, e 2 "^, ... , e' !2n/ > <M  l) ~] T . Unfortunately, this decomposition is not possible without
absolute knowledge of the noise. The p signal vectors e, contain the frequency infor
mation and are related to the eigenvectors of R„ bv v, = — ?= e, for i < p. The
eigendecomposition of R„ is
P M
= \ (). t + o 2 v )yjy? + > ajvfvf
(20)
The principal eigenvectors of R xx are a function of both the signal and noise. If no signal
is present the autocorrelation matrix is simply a diagonal matrix with the eigenvalues
equal to the variance of the noise [Ref. 1: pp. 422423]:
(21)
2
°v
The data generated from taking M samples of/? sinusoids in white noise can be
used to generate an autocorrelation matrix with the following properties:
• The theoretical autocorrelation matrix will be composed of a signal autocorrelation
matrix and a noise autocorrelation matrix.
• The signal autocorrelation matrix is not full rank if M> p.
• The p principal eigenvectors of the signal autocorrelation matrix may be used to
find the sinusoidal frequencies.
• The p principal eigenvectors of the signal autocorrelation matrix are identical to the
p principal eigenvectors of the total autocorrelation matrix.
The matrix will have a minimum eigenvalue ). = o]. [Ref. 1: pp. 422434]
Furthermore, the noise subspace eigenvectors are orthogonal to the signal
eigenvectors, or any linear combination of signal eigenvectors. For the theoretical
autocorrelation matrix, R M , all eigenvalues not associated with the signals are associated
with the noise and are identical in magnitude at ). = a] , the minimum eigenvalue of
R u . Thus,
R A/ V M = }mm y M (22)
The zeros of this minimum eigenvector polynomial
M\
I
y p+ ,(J+l)z J (23)
will have p zeros on the unit circle corresponding to signal frequencies [Ref. 4: pp.
335337, Ref. 5: pp. 371372].
For the simple case of M  1 complex sinusoids, the autocorrelation matrix R M
will have a single noise subspace eigenvector \ M with its associated eigenvalue / = a* ,
the minimum eigenvalue of R w . Thus, the resulting zeros correspond exactly to the
sinusoidal frequencies.
2. MUSIC
Multiple Signal Classification (MUSIC) is a form of a noise subspace fre
quency estimator, based on noise subspace eigenvectors with uniform weighting. The
MUSIC algorithm finds a pseudo spectral estimate from [Ref. 5: p. 373]
Puusictf) = JJT (24)
where e, = [1, e' 2 <, ... , e'< (M  1) ~] T . The advantage of this method is in its generalized
nature. There is no longer is a requirement for uniform spatial sampling. Any array
geometry will provide a solution. A welldesigned array will eliminate bearing ambigui
ties and provide unique solutions [Ref. 2: pp. 1928].
R.O. Schmidt [Refs. 2, 12. 13] has shown that a group of sensors excited by a
stationary point source emitter of known frequency will have a magnitude and phase
relationship (or correspondence) based on the DOA of the plane wave. This corre
spondence depends on the array geometry, as well as individual sensors characteristics
(i.e., directivity, gain, and frequency response), and may not be unique for that one di
rection of arrival.
By examining the signal in terms of an M dimensional complex field, where each
sensor provides an orthonormal axis, one can see that a single signal will result in one
steering vector. This steering vector describes the relationship between the individual
sensors in terms of phase and magnitude differences, thus for any unique signal fre
quency and direction of arrival there is one unique steering vector. The magnitude of
the vector will vary with time, but its direction in M space is a constant determined by
the amplitude and phase relationship of the sensors for that particular signal as illus
trated in Figure A.
The theory of superposition applies, thus with two signals present the instanta
neous received steering vector is a linear combination of the individual steering vectors.
The time varying steering vector will move in a plane that is spanned by the individual
steering vectors. Figure 5 shows the subspace plane spanned by two steering vectors.
The same idea can be expanded to a case of/? independent signals causing the steering
vector to move through a p dimensional subspace (as long as M> p).
Unfortunately, the steering vector may not determine the actual DOA uniquely.
In the example of a onedimensional linear array, a signal gi\es an infinite number of
Sensor 3
Sensor 1
x(0
Sensor 2
Figure 4. Steering vector for 3 sensors and 1 signal
Sensor 2
Sensor 1
Figure 5. Signal Subspace for 2 signals
DOAs that lie in a cone that is formed by revolving the actual DOA about the axis of
the array. Thus the array design and its manifold (expected response) will play a large
part in achieving optimum results. The array manifold describes the predicted response
of the array to any possible signal or combination of signals. The manifold may be es
timated analytically or through physical calibration.
Analytically describing an array is a complex mathematical procedure for all but
the simplest arrays (i.e., equally spaced sensors in a linear array or a 3 element
orthogonal array). It also assumes that absolute knowledge of the array geometry is
available  an assumption which can lead to errors if the differences are too large.
Calibration with test sources requires a known, low noise environment while the
test sources of each desired frequency are placed in a finite number of possible locations.
The estimated response to actual signals is an interpolation of these responses. Each set
of calibration parameters requires storage in memory; this results in overall massive
storage requirements for a comprehensive grid.
Several parameters such as the number of signals present, the directions of ar
rival of those signals, and the signal strengths can be determined from the signal sub
space information. More complex models, however, can be developed that can
determine signal frequency and polarization as described in References 2 and 13.
We will now describe the basic steps in the MUSIC algorithm for the DOA
problem. First, we sample the signals and compute the autocorrelation matrix of the
data R. Then, we determine the ordered set of eigenvalues and eigenvectors of R. In
particular, the eigenvectors associated with the minimum eigenvalues must be accurate.
In the theoretical case the signal eigenvalues are composed of signal strength (P, ) and
noise variance (o*) and the nonprincipal eigenvalues will all be identically a 2 v .
We now determine the number of signals by eigenvalue comparison. A simple
counting of the eigenvalues greater than a] will give the number of signals present in the
theoretical case. However, the sample autocorrelation estimates does not lead to such
simple results. Small power level signals may have small eigenvalues (perhaps smaller
than o] ). and the noise eigenvalues will not actually be identical but will group or cluster
near an approximation of a; . Methods of solution include likelihood ratio tests where
the gaps between the eigenvalues determine threshold placement [Ref. 2: pp. 7779].
Once we find the number of signals present we can determine the signal sub
space by the span of the first p eigenvectors
VkIv! v 2 ... v,] (25)
The signal nullspace is its orthogonal complement
y K c=[y p+] Vp+2 ... v M ] (26)
We now find the intersection of the signal subspace with the array manifold. The
intersection is given in equation 24 for the case of the uniformly spaced linear array. In
actuality the intersection is estimated with a "bestfit" method used to determine the
optimum result. In the nonlinear case the array manifold is much more difficult to rep
resent.
Two major areas of difficulty with the MUSIC algorithm are the calculation and
storage of the array manifold and performing the eigendecomposition of the
autocorrelation matrix that results from very large arrays.
3. ESPRIT
In Reference 14, Paulraj, Roy, and Kailath describe Estimation of Signal
Parameters via Rotational Invariance Techniques (ESPRIT), a subspace method which
utilizes two identical subarrays X and Y. A known distance called a displacement vector
separates the two parallel subarrays, but no rotation can be present. Each sensor in a
matching pair (doublet) must be identical, but knowledge of individual sensor and array
response characteristics is not required.
The N elements of both arrays are sampled simultaneously with the signal at
each sensor being given by
P
Sk(<) a i{ Q k) + n m
(27)
where the sampled signal at each sensor in a doublet differs only by a phase term and
additive noise. With /> signals present, s k is the wavefront at the reference sensor in the
X array, B k is the DOA relative to the displacement vector. a{6 k ) is the response of the
/th sensor in a subarray relative to the reference sensor in that array for a signal from
bearing 6 k , d is the magnitude of the displacement vector and n the noise term. In vector
notation the signals at the subarravs are
x(f) = Hs(0 + n x U)
\U) = HOs(r) + n v (/)
x(r) = lx 1 {t),x 2 {t),...,x m (t)j r t
y(0 = Oi(0. yiif),  , y m (t)l T ,
n x {!) = Zn Xl (t),n X2 (t),...,n Xm (t)l T ,
and n y (t) = [ w> . (r), n >2 (t), ... , n,^ 1
(28)
(29)
The vector of wavefronts at the reference sensor in array X is represented by s(/). All
displacements and phase differences are based on this sensor. The p by p diagonal ma
trix, O, contains the phase delays that occur with each set of doublets
<D = diag[^', ^...., J**] (30)
where 4> k = ^\~ sin 8 k , as shown in equation 13.
If R,, p is the signal autocorrelation matrix, the subarray autocorrelation matrix
is given by
R xx = UR pp U T + a 2 v l (31)
The cross correlation between subarrays is
R xy = HR„O r H r (32)
where H is the direction matrix whose columns are the direction vectors for the p
wavefronts
H0 k ) = UW, W h m (d k )Y (33)
After some manipulation [Ref. 15: pp. 251253], we obtain the matrix
r = (R„  W)  J**,
= HR pp H T  yHR pp ® T H T (34)
= UR pp {\  y<D r )H J
In general, there will be p eigenvalues of this matrix. But when y = e ,2 if sir,fl <, the /th row
of (I — y^ 7 ) becomes zero, leaving p\ eigenvalues. Each value of y where this occurs is
called a generalized eigenvalue (GE). Now, the DOAs can be determined by
d k = arcsm(^) (35)
Due to estimation errors in the calculation of the correlation matrices, some
error will be induced. Generally, the GE's will be moved off the unit circle and out from
the origin, but in the case of strong signals, they will still be easily discernible. It should
be noted that poor array design may result in possible ambiguities (similar to MUSIC).
Advantages to note over the MUSIC algorithm come with respect to the array
manifold. With ESPRIT, no manifold is required. The considerable calibration efforts
and storage requirements are nonexistent. However, the two subarrays must be identical
in all respects and must be positioned exactly parallel to each other [Ref. 14].
4. Other Subspace Methods
Large variance effects may be seen in the above methods due to the poor reli
ability of the estimation of the eigenvectors associated with the minimum eigenvalues
of the estimated autocorrelation matrix R„. A different method of spectral estimation
is the use of the principal components where only the eigenvectors associated with the
largest eigenvalues are used. Some methods have tried to minimize the effect of noise
by using R„ — o]\ but the estimation errors in both R„ and o 2 v have limited their suc
cess [Ref. 1: pp. 425]. Other spatial spectral methods include the parametric methods
such as AR and ML [Ref. 1: pp. 426427].
The structured matrix approximation of Kumaresan and Shaw [Ref. 16] uses K
snapshots in time of an N element uniformly spaced linear array with each snapshot in
time forming a column of a data matrix X, which is then approximated by X M in the least
squares sense. The bearings information is then calculated using the properties of the
Vandermonde matrix.
A combination of frequency domain beamforming and autoregressive modeling
techniques have been employed in Reference 17 to estimate the direction of arrival in a
multiple source localization situation using a planar array.
Halpeny and Childers [Ref. 18] use frequency wavenumber filters to break the
multiple wave case down to a succession of single wave problems.
Reference 19 explains an algorithm that uses nonnoise eigenvectors from a
covariance matrix to obtain high resolution direction of arrival for narrow band sources.
An adaptive beamforming method similar to a minimum energy approach is
decribed in Reference 20. The eigenstructure of the correlation matrix is analyzed and
the computations are modified to optimize resolution but at a cost of array gain.
III. THE LANCZOS ALGORITHM
Lanczos algorithms provide a method to find some of the extreme eigenvalues
(smallest or largest) and their associated eigenvectors from large matrices with fewer
operations than required in an entire matrix decomposition. The procedure is a
tridiagonalization of the original matrix based on an iteration developed by Cornelius
Lanczos [Ref. 21: pp. 49206]. Once the matrix is in a tridiagonalized form the
eigenvalues can be easily computed using a symmetric QR algorithm [Ref. 15: pp.
278279] or bisection (with or without Sturm sequencing) [Ref. 15: pp. 305308]. The
algorithm takes advantage of "minimized iterations" providing quick convergence to the
final tridiagonal matrix and avoiding accumulation of roundoff error. [Ref. 22]
The method fell from favor as a tridiagonalization technique with the advent of the
Givens and Householder transformations later on in the 1950s. Givens transformations
[Ref. 15: pp. 3847] use plane rotations (orthogonal matrices) to zero out undesired en
tries in the matrix undergoing tridiagonalization. The Householder methods [Ref. 15:
pp. 4347] use elementary reflectors to accomplish the same end. Both methods are in
herently stable, with the Householder method being slightly superior in both speed and
accuracy. Both methods outperformed Lanczos as a complete tridiagonalization proce
dure in the general case where all eigenvalues are required [Ref. 23: pp. 4243].
The real power of the Lanczos method however lies in obtaining the extreme values
quickly. The entire matrix need not be tridiagonalized before solutions start to converge.
Also, if the original matrix is sparse, the Lanczos method maintains that property, sav
ing even more computations. Thus for large matrices (order > 100) the Lanczos method
will converge on extreme eigenvalues in many fewer operations. Recently Tufts and
Melissinos [Ref. 24: pp. 4347] have derived a variation of the Lanczos method for high
resolution spectral estimation and showed that their method outperforms both the sin
gular value decomposition and the power method of principal eigenvector and
eigenvalue computation. Later in this chapter, it will be shown that storage require
ments are also minimized.
This chapter starts with an explanation of the classical Lanczos algorithm as devel
oped by Lanczos and modified by Paige [Ref. 25]. Then we will discuss the unorthodox
Lanczos algorithm of Cullum and Willoughby where no reorthogonalization is per
formed [Ref. 23]. Finally, the algorithm used in the direction of arrival problem are
discussed in detail. Also, we present some results of the algorithm in the form of spectral
estimation of multiple tones.
A. CLASSICAL LANCZOS AND ITS VARIANTS
The original algorithm was designed as a means to directly tridiagonalize the sym
metric matrix A. The development of Lanczos algorithm requires the knowledge of
Krylov sequences and subspaces. For an n by n matrix A and any arbitrary nonnull
n by 1 vector \\ the Krylov sequence of vectors is defined as:
v /+] = Avj = A'vj for k = 1, 2, ... , n (36)
In the above sequence there will be a vector, say v ro+1 , which can be expressed as a
linear combination of all the preceding vectors. The Krylov matrix of rank m is then
given by
V m = [v, v 2 v 3 ... vj = [v l5 Av,, A 2 v„..., A^vJ (37)
and the Krylov subspace is the space that spans the vectors [v„ \\. ... , vj,
A'" ; (A,v,) = span{v ls Av,, A 2 v ls ... , A^'vJ (38)
The columns of the n by m Krylov matrix V m are orthogonal. The tridiagonalization
of the given matrix A is then obtained as
T = v£AV OT , (39)
where T is an m by m tridiagonal matrix. Thus, the given matrix A can be
tridiagonalized provided we have an efficient way to compute the orthogonal matrix V ; .,
or to compute the elements of matrix T by performing the decomposition of equation
39 directly as proposed by Lanczos [Ref. 22, Ref. 15].
The most direct way of performing the tridiagonalization assumes that T = V 7 AV
where V = [v, v 2 ... v m ] . Note that A is a symmetric n by n matrix and V is an
n by m orthogonal matrix constructed from the Krylov sequence of vectors. The basic
recursion for a real n by n matrix A starts with a randomly generated unit Lanczos vec
tor v,. Define scalar /), = and v = 0. Define Lanczos vectors v, and elements tx, and
&_, for i = 1,2,..., M,
/^^Av.a^ftv, (40)
a, = v/Av,
(41)
(42)
This is referred to as the basic Lanczos iteration or recursion. The actual Lanczos ma
trix T, j = 1,2, ... , M, is a symmetric tridiagonal matrix with a, , 1 < /' <j, on the di
agonal, and /?,.! , 1 </<(/'— 1), above and below the diagonal.
(43)
«1
h
h
a 2
h
ft
•
•
.
o /?,
Each new Lanczos vector v /4 ., is obtained from orthogonalizing the vector Av, with
v, and v,_, . The elements a, and /?,., are the scalar coefficients that make up the Lanczos
matrix. The basic Lanczos recursion may be condensed into matrix notation by defining
^' m = ( v n v 2> ••• . V J Then from equations 40, 41, and 42
AV„
V„,T m
(44)
where e,„ is the coordinate vector whose mih component is 1 while the rest of the com
ponents are 0, J m is the m by m Lanczos matrix, the diagonal entries are
T m (k,k) = a k for \<k< m,
(45)
and the entries along the two subdiagonals on either side of the principal diagonal are
JJkJi + 1) = T m (k + \,k)  /J k+] for 1 < k < m  1 (46)
Note that A is never modified (unlike in Givens and Householder transformations), thus
advantage may be taken of any existing sparsity. Also, the only storage requirements
are for the three Lanczos vectors (n dimension), the Lanczos matrix T,, and the original
matrix A.
We summarize the actual steps:
• Use the basic recursion of equations 40, 41, and 42 to transform the symmetric
matrix A into a family of symmetric tridiagonal matrices T y , j= 1,2,...,.U.
• For m < M find the eigenvalues of T m , \x (also known as the Ritz values of TJ.
• Use some or all of these eigenvalues, n, as approximations of the eigenvalues of
A, / .
• For each eigenvalue ix Find a unit eigenvector u so that T m u = ^u .
The Ritz vector y is the approximation of the eigenvector of A. It is found from map
ping the eigenvector u of the Lanczos matrix.
y = V„u (47)
So the eigendecomposition of T m finally results in the Ritz pair (/x, y) which approxi
mates the eigenvalues and eigenvectors of A. [Ref. 23: pp. 3233., Ref. 15: pp. 322325.]
Unfortunately, the effects of finite precision arithmetic prevent the theoretical
Lanczos algorithm from working. The eigenvalues of A and the eigenvalues of T„ no
longer converge. Roundoff errors are partially to blame, but the dominant effect is from
the loss of orthogonality of the Lanczos vectors v, . From equation 40 it can be seen
that a small /?,_, will have great effect on v,.,. Paige showed that the algorithm runs
within allowable error constraints until it starts to converge on the first eigenvalue. At
this point />,.! becomes small and the Lanczos vectors lose orthogonality. The loss of
orthogonality is not random, however. It always occurs in the direction of the Ritz
vector y.
This trouble can be dealt with through reorthogonalization. But questions that we
must answer include:
• How much reorthogonalization is required?
• Where should it be performed?
• Reorthogonalize with respect to what?
Complete reorthogonalization is the first choice, inserting a Householder matrix
computation into the Lanczos algorithm. This enforces orthogonality among all the
Lanczos vectors and is effective at keeping the system stable. But the extra computa
tions it requires negate any advantage of the Lanczos algorithm, making the number of
computations on the same order as a complete decomposition. Numerous vectors have
to be stored and compared requiring many swaps into and out of storage. [Ref. 15: pp.
334335]
Selective reorthogonalization is clearly more efficient. The extra computations are
performed only if orthogonality checks indicate loss is imminent. Paige shows that the
loss of orthogonality occurs only when the algorithm converges on a Ritz pair. Thus,
instead of reorthogonalizing against every other Lanczos vector, using the less numerous
Ritz vectors will be sufficient. Storage is therefore minimized. Only when all eigenvalues
of the matrix are required does this method become computationally more expensive
than other techniques. [Ref. 26]
The last option is no reorthogonalization. The explanation above would seem to in
dicate that maintaining orthogonality is required. However by analyzing the causes and
effects of the original loss of orthogonalization one can insert corrections into the sol
ution to give valid eigenvalues and eigenvectors. This is the procedure that will be ex
amined in the next section.
One other property will be mentioned. Here the single vector Lanczos recursion
described above will not find duplicate eigenvalues of the matrix A. Parlett's proof uses
the power method to expand v to compute the columns of the Krylov matrix K"(v, A).
If J{\) is the smallest invariant subspace of R' 1 which contains v. then the projection of
A onto J{\) has simple eigenvalues. Also, the Krylov subspaces fill up J{\) and for
some / < /; we have
s P an{\] c K\k, v) c K\k, v) c ... a K\\, v) = K l+] (A, v) = J{\)
(48)
The result is that some eigenvectors of A may be missed, and any repeated eigenvalues
will not be detected. [Ref. 27: pp. 235239]
The multiple eigenvalue problem can be treated by using a Block Lanczos algorithm.
Instead of obtaining a tridiagonal matrix, the result is a banded block matrix, where the
diagonal is M„, an / by / matrix, and the above and below the principal diagonal are
upper triangular matrices B„ . Each block must be dimensioned / > p , where p is the
estimated multiplicity of the desired eigenvalue.
T,=
M]
Bi
B 2
M 2
Bf
B 3
•
•
B/
B
M
(49)
This is analogous to the general case [Ref. 15: pp. 337339]. The block Lanczos algo
rithm is briefly reviewed at the end of this chapter.
The above discussion assumes that the given matrix A is real and symmetric. Besides
the algorithms summarized in this section for a real and symmetric matrix case, there
are other general Lanczos algorithms proposed in the literature [Ref. 23] that are suitable
for Hermitian matrices, nonsymmetric matrices, and for rectangular matrices.
B. IMPLEMENTATION
The Lanczos phenomenon states that a few, not all, of the eigenvalues of a very
large matrix A can be computed using the Lanczos recursion with no
reorthogonalization. But to find most of the eigenvalues, the Lanczos matrix, T m , will
grow in size approaching that of A, causing the loss of orthogonality of the Lanczos
vectors. The loss of orthogonality results in many spurious eigenvalues, as well as extra
copies of good eigenvalues. In any case, a test is required to confirm either:
• a "found" eigenvalue is good, or
• the eigenvalue that appears is spurious.
Golub and Van Loan. Parlett, and Paige [Refs. 15, 27, and 28] describe procedures
that look at the eigenvalues for each T,„ as m is stepped up in size. All the eigenvalues
of l m are computed at each step. The good eigenvalues will repeat at each larger T m .
while the spurious eigenvalues jump around. If an eigenvalue does not match at con
secutive T„/s it may be considered spurious and thrown out. If a good eigenvalue is
mistakenly deleted (due to numerical roundoff), it can be counted on to reappear in the
next step.
Cullum and Willoughby [Ref. 23] take a different tack by developing a test to find
and eliminate bad eigenvalues and retain the rest. The advantage here is in utilizing the
machine's tolerance to drop bad eigenvalues, while not discarding good eigenvalues that
have yet to converge. As a result a larger spectrum is available sooner, even though it
may only be a rough estimate of where the eigenvalues will finally converge.
In practice, parts of the Lanczos recursion (equations 41 and 42) are replaced by
a, = v/lAv,./^,.,) (50)
and
fi w = IIAvjowftv,!! (51)
Computation of the element a, is a modified GramSchmidt orthogonalization proce
dure. The new /? /+ , is equivalent to the /?,^, of equation 42 but now it directly controls
the size of the Lanczos vector.
In what follows we describe two Lanczos algorithms, namely the single vector
Lanczos algorithm which is modified and analyzed by Paige [Ref. 28] and the block
Lanczos algorithm described by Cullum and Willoughby [Ref. 23]. Both algorithms have
been considered for the estimation of the directionsofarrival of multiple targets in noisy
environments in this thesis.
1. Single Vector Lanczos Algorithm
The first procedure to be described is the Paige's single vector Lanczos algo
rithm [Ref. 28] for real symmetric matrices. The single vector Lanczos procedure is one
of the most straightforward implementations of the theory. This procedure will find
some or many of the eigenvalues and eigenvectors of a real symmetric matrix A such that
Ax = /x. It will not detect repeated eigenvalues. However, it may be noted that for
many problems of interest in practice we do not have strictly multiple eigenvalues. For
example, in the directionofarrival problems the smallest eigenvalues of the
autocorrelation matrix corresponding to the noise associated with the target signals are
spread over a small range rather than coinciding on the same value [Ref. 2].
No reorthogonalization is performed as part of the single vector Lanczos al
gorithm. As mentioned earlier, the Lanczos vectors begin to lose their orthogonality
when we seek to estimate all or most of the eigenvalues of the real symmetric matrix A.
For the application under consideration, however, we are generally interested in only a
few of the minimum eigenvalues and the corresponding eigenvectors. It is mainly for this
reason that we have not attempted the complete or selective reorthogonalization of the
Lanczos vectors in this study.
Now we shall outline the basic steps involved in the single vector Lanczos al
gorithm. This is based on the recursion described by equations 40, 50 and 51. Based on
these equations Paige [Ref. 28] presented four different single vector algorithms. We
have adapted one of in this study. The complete eigenvalue/eigenvector problem actually
consists of three parts: (a) obtaining the tridiagonal matrix T m from the given symmetric
and real matrix A using Paige's recursion, (b) determining the smallest eigenvalues of
J m using the bisection method and Sturm sequencing, and (c) estimating the corre
sponding eigenvectors by computing the Ritz vectors. The details are presented in the
following.
Step 1: As shown in equation 43 the symmetric tridiagonal matrix T, has entries
a, and /?, along its principal, and the adjacent sub and super diagonals, respectively. The
following recursive expressions are then used to compute the entries of the tridiagonal
matrix, and also the Lanczos vectors v ; [Ref. 28]:
Initial conditions: v, is an arbitrary n by 1 vector such that v 2 = 1
u, = Av,
«o = fa =
for j = 1,2, ... , m
*j = v/u, (52)
W . = U;  Xffj
fij  li»>ll 2
V/ +]  YFjfij
u ;+] = Av y+1  fyvj
where w and u r are some intermediate vectors. The vector v, is obtained by filling its
entries with a random number sequence and then normalizing it with respect to its
Euclidean norm. Now T m is obtained by simply filling its entries as in equation 43. Note
that m 4 n in the above. One quick test to ensure that we have obtained a fairly accurate
estimate of T m is to compute the product v/v y . The result should be equal to 5 t] , where
S„ is the Kronecker delta function.
Step 2: The eigenvalues of the tridiagonal matrix T m , denoted by n p can be
computed using the bisection method and Sturm sequencing. Actually one could obtain
both eigenvalues and eigenvectors of T^ by employing such methods as the QR algo
rithm. However, when only a few eigenvalues are required, the bisection method seems
appropriate. For the given m by m matrix T m we define the characterstic polynomial
p m {n) = det (T m  al) (53)
which can be recursively computed as follows
; 7 o(^) = 1
p x {fi) = a x n
for j = 2, 3, ... , m
(54)
The roots of the polynomial p m {n) are the required eigenvalues. For our appli
cation, we are only interested in a small range of eigenvalues at the lowest end of the
eigenvalue spectrum. We make use of the Sturm sequencing property that the
eigenvalues of T,_, strictly separate those of T, [Ref. 15: pp. 305307] and implement the
following iteration:
a + b
b = v if Pm (a) Pm (b)<0 (55)
a = pl if p m (a)p m (b) >
and we repeat the above as long as  b — a \ > e(  b  +  a \ ), where c is the machine
unit roundoff error and \_a,b~] is the range of our required eigenvalues. Determining the
range of interest in our application may require some a priori knowledge about the signal
to noise ratios (SNR) and it may take a couple of iterations to do this. Some alternatives
to the iteration given in equation 55 are to use a straightforward polynomial root finder
and then pick the roots of interest, or to employ the well known LDU factorization,
both of which may not be computationally efficient.
Step 3: There are two ways to obtain the eigenvectors of A knowing its
eigenvalues, /,. Note that /i, are the estimates of A,. In the first method, we compute
the eigenvectors of T m , denoted by t , and then obtain the eigenvectors of A given by
x,. = Ljj (56)
where L m = fv, v 2 ... v m ] is the Lanczos matrix which ideally is the same as the
Krylov matrix of equation 37. Note that (li,, t) e T m .
The second method involves computing the Ritz vectors cither from the
Rayleigh quotient interation or by the orthogonal iteration. Here we assume that we
have good estimates of/, from the previous step, and proceed to obtain the eigenvector
Xj by minimizing the cost function
7= (A; 7 .I)x y  2 (57)
It can be shown that a simple minimization of J with respect to x, yields the Rayleigh
quotient of x.
r(xj) = Aj = ^4 (58)
XjXj
Therefore, given /, and using equation 58 we can formulate the Rayleigh quo
tient iteration as follows [Ref. 15, 27]
Initial condition: x is an arbitrary vector such that x  2 = 1
for k = 0, 1,2,...
r(*j) = ^L
tf* (59)
solve (A  r(x k )V)j k+] = x k for y k+]
**+] = h+ilhk+ih
where y kM is some intermediate vector. We stop the iteration when y A _, 2 converges to
a constant or when r(x k ) ^ ). k , one of the known eigenvalues. At each iteration step we
need to solve an nbx n system of equations in this method. One advantage with this
method, however, is that it converges very quickly. Besides the above iteration, some
other methods are outlined in References 15 , 23 and 27. We remark that if only a few
eigenvalues and eigenvectors (say, five) are required, it may be more direct to use
equation 56.
We now present an example of the ability of the single vector Lanczos algorithm
to estimate the direction of arrival, or to find spectral lines in noise, and the advantages
in extracting more than one eigenvalue and eigenvector in this process. We consider a
signal with three sinusoids present in noise
3
An) = y A t
cos(27ry>0 + n{n) (60)
where A, are the amplitudes of sinusoids, / are the normalized spatial frequencies
(0 <f : < 0.5 corresponding to .0 <6 <^) , and n(n) is the zero mean white noise with
a variance of a].
We have computed a 25 by 25 autocorrelation matrix of x(n), R„ , by using 100
data samples. We have used the covariance method for this purpose, hence R„ is real
and symmetric. The eigenvectors x, of R„ corresponding to the lowest eigenvalues are
computed by employing the single vector Lanczos algorithm. The power spectral density
estimates are computed as follows:
s%M
I
(61)
me ™
where x Jt are the elements of the jth eigenvector, x y .
Figure 6 and Figure 7 show the power spectral density (PSD) estimates of
equation 60 with an SNR of 10 dB for J = 1 and 2, respectively.
!
io'
1
10 2
IO" 3
u
A n /
1 !
io<
U
WW VIA ;
Figure 6. PSD of first eigenvector
Note that the index j indicates the increasing magnitude of the eigenvalues. Thus,
(/l,,x,) are the lowest eigenvalue and the corresponding eigenvector. In both figures we
have the peaks at the correct locations (9°, 27°, 63°). However, they both have spurious
smaller peaks at different locations. We can observe the same trend for the first five
eigenvectors as shown in Figure 8, where the spectral estimates are over laid on each
other.
Based on the above results one feels that we could employ some kind of aver
aging to get rid of the spurious peaks and improve the estimation performance. We have
implemented two such methods: the eigenvector averaging and the spectral averaging.
Figure 9 shows the result of the algebraic averaging of the first 5 eigenvectors, and
Figure 10 shows the result of the algebraic averaging of the spectral estimates of the
same eigenvectors. As seen from Figures 9 and 10, eigenvector averaging yields better
results than spectral averaging.
Figure 7. PSD of second eigenvector
Improved results, however, were obtained by using what is called spectral
multiplication which is obtained by taking the product of the individual spectra, given
by
r>
s«W = J J«(/)
(62)
where 7 is a predetermined number (7 <m< n). Figures 11 and 12 show the results of
spectral multiplication for 7 = 2 and 7=5, respectively. As can be seen in these fig
ures, using more spectra in equation 62 greatly improves the performance. Also, even
for 7=2, spectral multiplication outperforms the eigenvector averaging method.
Figure 8. Overlayed PSDs of first 5 eigenvector
In the remainder of the thesis, we have used spectral multiplication in prefer
ence to the eigenvector or spectral averaging. Figure 13 shows the multiplication of 5
spectra for the case when the SNR = dB. We notice a spurious peak around
6 = 74°. iMore spurious peaks are observed when the SNR is decreased to 5 dB (see
Figure 14) and Figure 15 shows the spectrum for the SNR but we have used the
eigenvectors 610 in this case. Improved performance is obtained as shown in
Figure 16 (J = 10) and Figure 17 (J = 15) by using more and more eigenvectors in the
spectral multiplication.
In all the above cases we always observed the signal spectral peaks at the right
places. The spurious peaks, however, did not appear at the same location as we used a
different eigenvector to compute the spectrum, S^(/).
Figure 9. Eigenvector averaging
Figure 10. Spectral averaging
10"
io»
10z
■f ,0 ~ 3
c
a 10" 4
a
oo 10 «
J
^A
u
w
108
"V
vy
10»
\J
20
40
Ben
60
rings
80 1
30
Figure 11. Spectral product for 2 eigenvectors
10°
10s
c
o
3 10 ~ 12
V
1 
A
a.
io'»
u
U
M/
w
10" 24
c
20
40 60
80
1(
)0
Bearings
Figure 12. Spectral product for 5 eigenvectors
36
io'
102
■f ,0 ~ 3
c
a 10"<
H IO" 5
IO" 6
io°
l
lo io
C
,
II
20 40
60
80 100
De
rings
Figure 13. Spectral product for 5 eigenvectors, dB
10'
102
5 103
c
Q IO" 4
tl ,0 " S
W 10«
O IO" 7
io«
io»
Ill
io'°
1 .A 1 1 /. 11 1.
3 20 40 60 80 1
30
Bearings
Figure 14. Spectral product for 5 eigenvectors, 5 dB: Using second through sixth
eigenvectors
 10« "
'
io'
io 2
'
2 IO" 3
a
£ io<
II
% io 5
V
OT 10«
1
o IO" 7
ft
io°
I
io»
 1

io'°
1
,
1 .
J
Ill,
20 40
60
80 100
Bear
MRS
Figure 15. Spectral product for 5 eigenvectors, 5 (IB: Using sixth through tenth
eigenvectors
10"
io'
IO" 2
§ IO" 3
Q 10"<
x> l0 " s
W 108
5 ,o
IO" 8
IO" 9
1
1010
20
40
Bear
60
ings
80
100
Figure 16. Spectral product for 10 eigenvectors, 5 dB
10"
io'

io«

•f 10 ~ 3
c
a 10<

£ 10»
W , 6
o io'
IO" 8
10»
IO"'"
(
I 20 40 60
80 100
Bearings
Figure 17. Spectral product for 15 eigenvectors, 5 dB
2. Other Methods
The single vector Lanczos algorithm will not determine that repeating
eigenvalues exist, thus it cannot find the corresponding eigenvectors. The subspace that
results has an incomplete basis as it is described only by the eigenvectors that are com
puted. The solution is to use a block method that is analogous to the single vector
Lanczos algorithm. As we mentioned earlier, the block form of the Lanczos algorithm
does find eigenvectors with multiplicity p as long as the blocks are dimensioned l>p.
We attempted to incorporate the Cullum and Willoughby hybrid block Lanczos algo
rithm (Ref. 29 Chapter 8) into our direction of arrival model. We postulated that it
would be desirable to compute a few of the extreme smallest repeating eigenvalues and
their respective eigenvectors. However we were never able to get the program to reliably
compute good eigenvalues and eigenvectors for the autocorrelation matrix. This has not
posed a problem for our model as the covariance matrix does not appear to have re
peating eigenvalues, but a larger order matrix may indeed include duplicating noise
eigenvalues and require an algorithm that will accurately operate with that perturbation.
The algorithm we attempted to use is actually a hybrid approach to finding the
eigenvalues and eigenvectors of a real symmetric matrix A. For insight into the problem
look at the block analogy of equations 40, 41, and 42. Define matrices
B t = and V„ = 0. The « by q matrix \ x has columns that are orthonormal random
vectors. The value of q must be greater than or equal to the number of eigenvalues to
be found.
for / = 2, 3, ... , 5
t (63)
M, = V/CAVjV^B/) (64)
V W B W  P, (65)
The matrix B,,, is a modified GramSchmidt orthonormalization of the columns of P,,
Also note that the matrix M, correspond to the pc's of the single vector Lanczos.
The block analogy to the Krylov subspace approach can be performed with
/^(A.V,) = span{\ lt AV lt A 2 V,, ... , A s ~ ] \] (66)
The blocks V„ for j = 1, 2, ... , 5 form the orthonormal basis of the Krylov subspace.
It can be shown that for a symmetric n by n matrix A and an orthonormal
n by q starting matrix V,, that the block recursion equations 63, 64, and 65 will generate
blocks V 2) V 3 , ... , \ s where qs < n. It is these blocks that form an orthonormal basis
of the subspace A''(A, V,) . In much the same way as the single vector Lanczos algorithm
generates the tridiagonal Lanczos matrix, the block variant generates blocks, but these
are now nontridiagonal. At the end of each iteration the Lanczos matrix is of the form
A,
«
Mj AM
b 2
«2
&
M r AMj
h
<*3
04
0C 4
ft
(67)
The submatrix MfAM consists of the reorthogonalized terms and M, is the portion of
the first block that is not generating descendants. Ritz vectors are computed on every
iteration and are used as the starting blocks for the next iteration. Each block is required
to be reorthogonalized with respect to the all the vectors in first block which is not being
allowed to generate descendants. It is apparent that the block procedure requires a great
amount of storage and is very computationally intensive.
IV. RESULTS
Using the Lanczos algorithm it is possible to find some of the eigenvalues and
eigenvectors of a matrix without going through an entire matrix decomposition. The
smallest eigenvalue of the autocorrelation matrix and its corresponding eigenvector will
have the required spectral information to determine a source's bearing (direction of ar
rival) from an array. Multiplying several of the resultant eigenvectors' power spectral
densities will tend to reinforce the true spectral peak and zero out spurious peaks that
do not occur with every eigenvector.
The problem with finding the split between the noise and signal eigenvalues disap
pears as only a few of the smallest eigenvalues of a large matrix (in relation to the
number of sources) are used.
A. APPROACH
The received signal is modeled by sinusoids at normalized spatial frequencies pro
portional to their bearings from endfire (.0 = 0° , .5 = 90° ). The sum of these
sinusoids is sampled at a rate based on the interelement spacing of A m J2. Thus a source
at endfire is sampled at the Xyquist rate and the sample rate increases as the bearing
shifts toward the arrav broadside. The simulation uses
T
1
ss(n) = ) A cos(2tt/«) + n{n) (68)
where ss{n) is the instaneous excitation for the sensor at location n, A is the amplitude
of each of the T signals, f, is the normalized spatial frequency of the fth source (de
pendent upon bearing), and n{n) is white gaussian noise. The relationship between A
and the noise variance a] is determined by the desired signal to noise ratio (SNR), where
SNR = lOlogf y" J (69)
The experiment consists of simulating a linear array with equally spaced sensors re
ceiving signals of known temporal frequency from various bearings. One possible
physical implementation would place a bank of bandpass filters on each sensor with the
outputs from each similar filter tied into a correlator. Advantages of this method include
the processing gain found by prefiltering the noise and simple parallel implementation
with separate channels for each frequency band. The lowering of the noise bandwidth
will raise the SXR at the correlator. As more filters are used (smaller bandwidth) the
noise power decreases and the SXR is increased. The algorithm creates an
autocorrelation matrix with the output of the correlator. The Lanczos tridiagonalization
and eigendecomposition provides the eigenvectors that are estimates of the spatial PSD.
The PSD corresponds to the sources directions of arrival. A possible implementation is
shown in Figure 18.
B. EXPERIMENT SET UP
The first three cases show the effect of different signal strengths on the ability to
accurately determine the number of targets and the bearing resolution for various di
rections and target spacing. In each of these cases, the number of sensors is 100, a ma
trix size of 25 is used and 15 iterations (the number of eignevalues/ eigenvectors found)
are performed. Case 4 uses three 5 dB sources at 18°, 36° and 41.4° to illustrate the ef
fects on changing the number of sensors (samples), the size of the autocorrelation ma
trix, and the number of eigenvectors used. The noise is randomly generated white
gaussian noise with a standard deviation selected to provide the desired SXR. Each
figure shows, (a) the PSD of selected eigenvectors overlayed and plotted versus bearing
and (b) the product of selected PSDs of those eigenvectors.
Case 1 is with all sources at a signal strength oC 5 dB. Figure 19 shows results from
the second through sixth eigenvectors of a single source at 1S°. Xote that some of the
eigenvectors have individual peaks as high as the true signal peak, but only at the true
bearing do all have a common peak. Figure 20 illustrates the other end of the spectrum,
at 81°. Once again the second through sixth eigenvectors are overlayed to show that the
correct bearing is consistently displayed, but in this case one eigenvector has an indi
vidual peak higher than the true signal peak. The product of these PSDs provides suf
ficient resolution. Figure 21 has two closely spaced sources at 36° and 38°. Resolution
is achieved but the PSD product of the second through sixth eigenvectors shows a spu
rious peak near 75°. Figure 22 is from three sources at 0°, 36° and 88.2°. The individual
eigenvector PSDs clearly show the excellent performance at broadside.
Case 2 lowers the signal strength of all sources to dB. Figure 23 shows results
from the second through sixth eigenvectors of a single source at 18°. Many more peaks
I
I
L
SENSORS
BANDPASS
FILTERS
CORRELATORS
E I GE NDECOMPOS I T I ON
SPECTRAL
ESTIMATOR
bearing:
Figure 18. A physical implementation
are visible in the PSD product, making the decision of how many targets more difficult,
figure 24 shows that at the the other end of the spectrum (at 81°) the situation is
20 40 60 80
Bearings
Figure 19. Case 1 5 dB, 1 target at 18 °: (a) overlay of PSDs from second
through sixth eigenvectors, (b) product of the above PSDs
10"
io'

102

£ io 3
c
D IO" 4
" io«
V
10" 8
1
"
10»
"
II
i
1
40 60 80 100
Bearings
Figure 20. Case 1 5 (IB, 1 target at 81 °: (a) overlay of PSDs from second
through sixth eigenvectors, (b) product of the above PSDs
Iff
AA./
MMju. 11 '
Bearings
Bearings
00 100
10'
1
r
102
£ 103

c
Q 1(H

b io5
.
c/> , e
I ,0
1
io°
1 1 1
io»

,

. W.ki
80 100
Figure 21. Case 1 5 (IB, 2 targets at 36 ° and 38 °: (a) overlay of PSDs from third
through seventh eigenvectors, (b) product of the above PSDs
00 too
n
Ueari
ngs
I0«
1 i
c
Q
1
\
J I"' 2
a
A ;
o
1
\ A ;
*" oin
w v
VlT !
20 40 GO 00 J 00
Bearings
Figure 22. Case 1 5 dB, 3 targets at °, 36 ° and 88.2 °: (a) overlay or PSDs
from second through sixth eigenvectors, (b) product of the above PSDs
slightly worse (due to a lower sampling rate). The overlays of the second through sixth
eigenvectors show that the correct bearing is consistently displayed, but in this case
enough spurious peaks reinforce one another, resulting in the PSD product that has not
zeroed out the bad peaks. Figure 25 illustrates that the proper choice of eigenvectors
will resolve this problem. Here the PSDs of the first through fifth eigenvectors are used,
giving a product that is easier to determine correctly. Figure 26 shows the dB case for
two close spaced sources at 36° and 38°. Resolution is achieved but the PSD product
of the first through fifth eigenvectors shows several spurious peaks, including the same
one as in Case 1 near 75°. Figure 27 is the three source example at dB. The individual
eigenvector PSDs are repeating at the proper bearings but the performance at broadside
is resulting in the product at the other bearings actually being driven down.
A signal strength of 5 dB is used for Case 3. Figure 2S shows results from the first
10 eigenvectors of a single source at 18°. Using more good eigenvectors increases the
likelihood that all spurious peaks will be diminished. Figure 29 illustrates results at the
other end of the spectrum (at 81°). Resolution is looked at in Figure 30. At 5 dB the
algorithm cannot separate the two close spaced sources at 36° and 38°. A number of
spurious peaks are higher than the bump at 36° making it impossible to accurately de
termine the number of sources as well as both locations. Resolution is tried again in
Figure 31 with 2 sources at 36° and 40° . Using five eigenvector PSD products produces
good results. Figure 32 shows 3 sources at 5 dB. Good performance is seen both in
the overlays and in the PSD product.
Case 4 starts with 100 sensors and a 25 by 25 covariance matrix shown in
Figure 33. As the number of sensors is decreased and the number of eigenvectors used
is held constant, more spurious peaks start to occur (Figure 34 and Figure 35).
Figure 36 shows the spectral improvement as more eigenvectors are used. As the
number of sensors is decreased to 40 (Figure 37) and seven eigenvectors are used, the
results are still acceptable. Using only 30 sensors we can no longer see resolve between
the two closely spaced sources. With a sufficient number of eigenvectors no spurious
peaks are present but the true number of targets is nondeterminable (Figure 38 and
Figure 39).
10'
lO 2

1 ,0 ~ 3
c
o 10<
1)
a
W ()0
1 ,u ~ 7
ft
10" e
1
io°
H
inio
1
. 11.11. 11
20 40 60 00 100
Bearings
Figure 23. Case 2 dB, 1 target at 18 °: (a) overlay of PSDs from second
through sixth eigenvectors, (b) product of the above PSDs
102
103
Li
Bearings
Figure 24. Case 2 cIB, 1 target at 81 °: (a) overlay of PSDs from first through
fifth eigenvectors, (b) product of the second through sixth eigenvector
PSDs
10°
loJ


102


£ io3
to
c

a 10 4
Power Spectral
o o o

1

108
loe
•
V

lOio
, 1 1 , 1 ,
C
> 20 40 60 80 100
Bearings
Figure 25. Case 2 dB, 1 target at 81 °: product of the first through fifth
eigenvector PSDs
Figure 26. Case 2 dB, 2 targets at 36 ° and 38 °: (a) overlay of PSDs from first
through third eigenvectors, (b) product of the first through fifth
eigenvector PSDs
40 60 80 100
Bearings
Figure 27. Case 2 dB, 3 targets at °, 36 ° and 88.2 °: (a) overlay of PSDs
from second through sixth eigenvectors, (b) product of the above PSDs
i no
i. ,rr
!
,.,
Ul
.
c
,1/11
a.
W , 3
*
o
cu
 \ J V y V/
io"
\ \1 \
=
I0" s
xj v \
20
40 60 00
Hearings
100
10 n
io«
102
£ I0" 3
c
Q lO 1
Q.
w 10" 6
1 ..
io°
100
1

1010
1 1
. II .
20
40 60 80
Bearings
100
Figure 28. Case 3 5 dD, 1 target at 18 °: (a) overlay of PSDs from first through
sixth eigenvectors, (b) product of the first 10 eigenvector PSDs
io'

Id" 2

£ io3
c
V
Q IO" 4

5 IO" 5
D
a
« 108
O .07

IO"
io°
IOI0
1
Bearings
Figure 29. Case 3 5 (IB, 1 target at 81 °: (a) overlay of PSDs from first through
sixth eigenvectors, (b) product of the first 10 eigenvector PSDs
Bearings
10°
Bear
ings
io'
102
.
£ 10" 3
c
*>
Q IO" 4
t> l0 " S
a
w ]0 _e

k '■>'
J
io°
io°
II
\\l
] IO
I
1 111
1
Figure 30. Case 3 5 dB, 2 targets at 36 ° and 38 °: (a) overlay of PSDs from first
through sixth eigenvectors, (b) product of the above PSDs
10° r
[
io'
llll
I
c
Q \0 2
nVM/M
;
D.
to , a
t
o
a
i
io 1
i
IO" 5
c
20 40 60
00
100
Bearings
10°
io'
I
IO' 2

*:' io 3
.
a
J
w , o
ID
O IO"'
a.
io°
io'
,oio
(
) 20 40 60
00
1
50
Bearings
Figure 31. Case 3 5 dB, 2 targets at 36 ° and 40 °: (a) overlay of PSDs from first
through sixth eigenvectors, (b) product of the first 5 eigenvector PSDs
Bear
ings
io>
lO 2
.
£ 103
c
V
Q 10"<
1
a
m io«
u
I)
 10 " 7
"
108

io»
inio
1,
1
• 1
Bearings
Figure 32. Case 3 5 dB, 3 targets at °, 36 ° and 88.2 °: (a) overlay of PSDs
from first through sixth eigenvectors, (b) product of the first 10 PSDs
Figure 33. Case 4 5 dB, 3 targets at 18 °, 36 ° and 41.4 °: 100 sensors, (a) PSDs
of second through sixth eigenvectors, (b) product of the above PSDs
Figure 34. Case 4 5 dB, 3 targets at 18 °, 36 ° and 41.4 °: 75 sensors, (a) PSDs
of second through sixth eigenvectors, (b) product of the above PSDs
i
io'
:
£>
:
a
a 10" 2
«" 103
t
w / « n fl A 
CL
1 J \ j \J \V a I M /I I\
io<
J 1 \ 1 \ VahI 1 (An
in5
vy \ / Vy Y I I
Searings
io'
102

c
O 10"<
.
i io»
.
D.
W l0 8
1 10 " 7

IO"
V

io»
lOio
Lkl
.
Figure 35. Case 4 5 dB, 3 targets at 18 °, 36 ° and 41.4 °: 50 sensors, (a) PSDs
of second through sixth eigenvectors, (b) product of the above PSDs
62
10°
ioi

lO 2

£ 103
w
c
<v
Q 10" 4
;
J
To
tl 105
W io6
a>
O 107
;
10~ 8
•
J
_
109
•
.
lOiol

1
20 40
60 80 100
E
learings
Figure 36. Case 4 5 riB, 3 targets at 18 °, 36 ° and 41.4 °: Product of the fust
through eighth eigenvector PSDs
20 40 60 80 100
Bearing?
Figure 37. Case 4 5 (IB, 3 targets at 18 °, 36 ° and 41.4 °: 40 sensors, 20 by 20
matrix, (a)PSDs of first through seventh eigenvectors, (b) product of
the above PSDs
64
10"
10'
c
Q 102
1
:
1 10 " 3
t
a.

io<
\
io 5
(
) 20
40 60
80
100
10«
io'
Bearings
IO" 2
t 10 ~ 3
a
a io«
1 io»
00 io°
IO" 8
10»
10iu
J 20
40 60
00
1
JO
Bearings
Figure 38. Case 4 5 dB, 3 targets at 18 °, 36 ° and 41.4 °: 30 sensors, 20 by 20
matrix, (a)PSDs of first through fifth eigenvectors, (b)Product of the
first through fourth eigenvector PSDs
£ io' :
20 40 60 80 100
Bearings
60 80
Figure 39. Case 4 5 dB, 3 targets at 18 °, 36 ° and 41.4 °: 30 sensors, 20 by 20
matrix, (a)Product of the first through fifth eigenvector PSDs,
(b)Product of the first through twelfth eigenvector PSDs
V. CONCLUSIONS AND RECOMMENDATIONS
The results plotted in Chapter IV indicate that the eigenvectors found using the
Lanczos algorithm are sufficiently accurate to determine the spectrum. Although no
direct comparisons with other eigendecomposition methods are performed, the theory
indicates that many fewer operations are required. We handle the other difficulty of
conventional subspace methods by using only a few of the eigenvectors associated with
the minimum eigenvalues of the autocorrelation matrix. No estimation of the noise
subspace dimension is required or performed.
This theory may be applied to any system requiring rapid decomposition of the
correlation matrix. Examples include phased array radar and passive acoustic arrays
[Refs. 30, 31]. Reference 32 details an experimental system using the MUSIC algorithm
for multiple source direction finding.
The following areas are recommended for future study.
• Use of the products of multiple spectra apparently resulted in good detection at low
SNR. More research in this area to determine a physical interpretation of this
method is required.
• Analysis and comparison of the results in terms of computational speed and accu
racy with other eigendecomposition methods should be performed to find the true
results.
• The Lanczos algorithm developed uses no reorthogonalization nor will it find re
peating eigenvalues. Other forms of the Lanczos algorithm are available. Com
parisons between these different methods to determine accuracy and speed may
lead to more optimum results.
• A more detailed model should be developed that will simulate an array with a bank
of bandpass filters to better forecast results of a physical implementation (as seen
in Figure IS of Chapter IV).
• A method which implements the algorithm in parallel fashion may be tried. Using
one long linear array, several overlapping subarrays may be used to simultaneously
create several autocorrelation matrices. The algorithm may be applied to these
matrices in parallel. It is predicted that the greater number of available
eigenvectors will more properly describe the noise subspace and therefore more
accurately estimate the spectrum.
LIST OF REFERENCES
1. Kay, Steven M., Modern Spectral Estimation, PrenticeHall, Englewood Cliffs, 1987
2. Schmidt, Ralph Otto, "A signal subspace approach to multiple emitter location and
spectral estimation," Ph.D. Thesis, Stanford University, Stanford, California, Nov.
1981
3. Kay, Steven M., and Stanley Lawrence Marple, Jr., "Spectrum analysisa modern
perspective," Proceedings of the IEEE, Vol. 69, no. 11, Nov. 1981, pp. 13801419
4. Orphanidis, Sophocles J., Optimum Signal Processing, An Introduction, 2d ed.
Macmillan, New York, 1985
5. Marple. S. Lawrence, Jr., Digital Spectral Analysis with Applications, PrenticeHall,
Englewood Cliffs, 1987
6. Kay, Steven, and Cedric Demeure, "The highresolution spectrum estimator—a
subjective entity," Proceedings of the IEEE, Vol. 72, no. 12, Dec. 1984, pp.
1S151816
7. Harris. Fredric J., "On the use of windows for harmonic analysis with the discrete
fourier transform," Proceedings of the IEEE, Vol. 66, no. 1. Jan. 1978, pp. 5183
8. Nuttall. Albert H., "Some windows with very good sidelobe behavior," IEEE Trans,
on ASSP, Vol. ASSP29, no.l, Feb. 1981, pp. 8491
9. Van Veen, Barry D., and Richard A. Roberts, "Partially adaptive beamformer de
sign via output power minimization," IEEE Trans, on ASSP, Vol. ASSP35, no. 11,
Nov. 19S7, pp. 15241532
10. Van Veen. Barry D., and Kevin M. Buckley, "Beamforming: a versatile approach
to spatial filtering," IEEE ASSP Magazine, Vol. 5, no. 2, Apr. 1988, pp. 424
11. Dudgeon, Dan E., and Russell M. Mersereau, Multidimensional Digital Signal
Processing, PrenticeHall, Englewood Cliffs, 1984
12. Schmidt, Ralph O., "New mathematical tools in direction finding and spectral
analysis," Proceedings of SPIE, Real Time Signal Processing VI, Vol. 431, Aug.
2325, 19S3, San Diego, CA, pp. 719
13. Schmidt, Ralph O., "Multiple emitter location and signal parameter estimation,"
IEEE Trans, on Ant. and Prop., Vol. AP34, no. 3, Mar. 1986, pp. 276280
14. Paulraj, A., R. Roy, and T. Kailath, "Estimation of Signal Parameters via Rota
tional Invariance TechniquesESPRIT," Proceedings of the J 9th Asilomar Conf. on
Ckts.. Syst., and Comp., Pacific Grove, 1986, pp. 8389
15. Golub, Gene 11., and Charles F. Van Loan, Matrix Computations, Johns Hopkins
Press. Baltimore, 19S3
16. Kumaresan. Ramdas, and Arnab K. Shaw. "Superresolution by structured matrix
approximation," IEEE Trans, on ASSP, Vol. ASSP36, no. 1, Jan. 1988, pp. 3444
17. Nikitakos. Nikitas V., "A comparison of two frequency domain adaptive
beamforming algorithms for sonar signal processing," Master's thesis, Naval Post
graduate School. Monterey, California, 1988
IS. Halpeny, Owen S.. and Donald G. Childers, "Composite wavefront decomposition
via multidimensional digital filtering of array data," IEEE Trans, on Ckis. and
Syst., Vol. CAS22, Jun. 1975, pp. 2762S6
19. Cadzow, James A., "A high resolution Direction ofArrival algorithm for narrow
band coherent and incoherent sources," IEEE Trans, on ASSP, Vol. 36, no. 7, Jul.
19S8, pp. 965979
20. Johnson, Don H., and Stuart R. DeGraaf, "Improving the resolution bearing in
passive sonar arrays by eigenvalue analysis," IEEE Trans, on ASSP, Vol. ASSP30,
no. 4. Aug. 1982. pp. 638647
21. Lanczos, Cornelius, Applied Analysis, PrenticeHall, Englewood Cliffs, 1956
22. Lanczos, Cornelius, "An iteration method for the solution of the eigenvalue prob
lem of linear differential and integral operators," Journal of Research of the National
Bureau of Standards, Vol. 45, no. 4, Oct. 1950, pp. 255282
23. Cullum, Jane K., and Willoughby, R.A., Lanczos Algorithms for Large Symmetric
Eigenvalue Computations, Vol. I, Birkhauser, Boston, 1985
24. Tufts, D. and Melissinos, CD., "Simple, effective computation of principal
eigenvectors and their eigenvalues and application to highresolution estimation of
frequencies," IEEE Trans, on ASSP, v. ASSP34, no. 5, Oct. 19S6, pp. 10461053
25. Paige, C. C, "Computational variants of the Lanczos method for the
eigenproblem," J. Inst. Math. Appl, 10, 1972. pp. 373381
26. Paige, C. C, "The computation of eigenvalues and eigenvectors of very large sparse
matrices," Ph.D. Thesis, University of London, London, England, 1971
27. Parlett. B. X., The Symmetric Eigenvalue Problem, PrenticeHall, Englewood Cliffs,
19 SO
28. Paige, C. C, "Accuracy and effectiveness of the Lanczos algorithm for the sym
metric eigenproblem," Linear Algebra and its Applications, 34, American Elsevier
Publishing (1980), pp. 235258
29. Cullum, Jane K., and Willoughby, R.A., Lanczos Algorithms for Large Symmetric
Eigenvalue Computations, Vol. 2, Birkhauser, Boston, 1985
30. Nickel, L\, "Angular superresolution with phased array radar: a review of algo
rithms and operational constraints," IEE Proc, Vol. 134, PT. F, No. 1, Feb. 1987,
pp. 5359
31. Malpass, Robert L., "ATAS: Big ASW help for small ships," U.S. Naval Institute
Proceedings, Vol. 1 15 1 1031, Jan. 1989, pp. 107108
32. Schmidt, Ralph O., and Raymond E. Franks, "Multiple source DF signal process
ing: an experimental system," IEEE Trans, on Am. and Prop., Vol. AP34, no. 3,
Mar. 1986. pp. 281290
INITIAL DISTRIBUTION LIST
No. Copies
1. Defense Technical Information Center 2
Cameron Station
Alexandria, VA 223046145
2. Library, Code 0142 2
Naval Postgraduate School
Monterey, CA 939435002
3. Chairman, Code 62 1
Department of Electrical and Computer Engineering
Naval Postgraduate School
Monterey. CA 939435000
4. Professor Murali Tummala, Code 62Tu 3
Department of Electrical and Computer Engineering
Naval Postgraduate School
Monterey, CA 939435000
5. Professor C. W. Therrien, Code 62Ti 1
Department of Electrical and Computer Engineering
Naval Postgraduate School
Monterey, CA 939435000
6. Professor Ralph Hippenstiel, Code 62Hi 1
Department of Electrical and Computer Engineering
Naval Postgraduate School
Monterey, CA 939435000
7. Jeffrey M. Speiser, Code 743 1
Naval" Ocean Svstems Center
San Diego, CA 92 1525000
S. Dr. Rabi Madan 1
Office of Naval Research
Code 1114
800 North Quincv Street
Arlington, VA 222175000
9. Daniel E. Gear 4
104 Mervine Street
Monterey, CA 93940
10. Yons Joo Kim 1
SMC 1017
Naval Postgraduate School
Monterey, CA 939435000
Thesis
G2558 Gear
DOA estimation by eigen
decomposition using single
vector Lanczos algorithm.
Gear
jOA estimation by eigen
2Compos ? t'on using single
vector Lanczos algorithm.