U.S. patent application number 09/947026 was filed with the patent office on 2003-05-08 for digital beamforming radar system and method with super-resolution multiple jammer location.
This patent application is currently assigned to Lockheed Martin Corporation. Invention is credited to Yu, Kai-Bor.
Application Number | 20030085832 09/947026 |
Document ID | / |
Family ID | 25485390 |
Filed Date | 2003-05-08 |
United States Patent
Application |
20030085832 |
Kind Code |
A1 |
Yu, Kai-Bor |
May 8, 2003 |
DIGITAL BEAMFORMING RADAR SYSTEM AND METHOD WITH SUPER-RESOLUTION
MULTIPLE JAMMER LOCATION
Abstract
In a radar system, sampled aperture data are received from an
antenna array. The sampled aperture data include data that do not
correspond to echo returns from a beam transmitted by the antenna.
A covariance matrix is generating using the sampled aperture data.
An eigenvalue decomposition is performed on the covariance matrix.
A direction of arrival is determined from which at least one jammer
is transmitting a signal included in the sampled aperture data,
based on the eigenvalue decomposition.
Inventors: |
Yu, Kai-Bor; (Niskayuna,
NY) |
Correspondence
Address: |
DUANE MORRIS LLP
100 COLLEGE ROAD WEST, SUITE 100
PRINCETON
NJ
08540-6604
US
|
Assignee: |
Lockheed Martin Corporation
|
Family ID: |
25485390 |
Appl. No.: |
09/947026 |
Filed: |
September 5, 2001 |
Current U.S.
Class: |
342/16 ; 342/157;
342/158; 342/17; 342/90; 342/96; 342/97 |
Current CPC
Class: |
G01S 3/74 20130101; G01S
7/36 20130101; G01S 13/68 20130101 |
Class at
Publication: |
342/16 ; 342/17;
342/90; 342/96; 342/97; 342/157; 342/158 |
International
Class: |
G01S 007/36 |
Claims
What is claimed is:
1. A method for locating a radar jammer, comprising the steps of:
(a) sampling aperture data received by an antenna array, wherein
the sampled aperture data include data that do not correspond to
echo returns from a beam transmitted by the antenna; (b) generating
a covariance matrix using the sampled aperture data; (c) performing
an eigenvalue decomposition on the covariance matrix; (d)
determining a direction of arrival from which at least one jammer
is transmitting a signal included in the sampled aperture data,
based on the eigenvalue decomposition.
2. The method of claim 1, further comprising the step of providing
the direction of arrival to a tracker.
3. The method of claim 1, wherein step (d) includes determining an
angular location of a jammer that is not within a mainlobe of the
antenna array.
4. The method of claim 1, wherein step (d) includes determining
angular locations of a plurality of jammers that are located within
one or more sidelobes of the antenna array.
5. The method of claim 4, wherein step (d) includes determining an
angular location of an additional jammer that is within a mainlobe
of the antenna array.
6. The method of claim 5, wherein step (d) includes determining
angular locations of a further plurality of jammers that are within
the mainlobe of the antenna array.
7. The method of claim 1, wherein step (d) includes determining
angular locations of a plurality of jammers that are within a
mainlobe of the antenna array.
8. The method of claim 1, wherein steps (a) through (d) are
performed by a radar system that is operating in surveillance
mode.
9. The method of claim 1, further comprising: (e) updating the
covariance matrix repeatedly based on additional sampled aperture
data; (f) updating the eigenvalue decomposition each time the
covariance matrix is updated; (g) determining a number of jammers
within a field of view of the antenna based on the updated
eigenvalue decomposition.
10. The method of claim 9, wherein the number of jammers within the
field of view is determined based on a number of significant
eigenvalues of the covariance matrix.
11. The method of claim 9, wherein the updates of step (e) are
performed at regular intervals.
12. The method of claim 11, wherein the regular intervals have a
period of between about one millisecond and about ten
milliseconds.
13. The method of claim 9, wherein a super-resolution subspace
algorithm is used to estimate a direction of arrival of signals
from each of the plurality of jammers.
14. The method of claim 9, wherein step (e) includes using a
sliding window to update the covariance matrix.
15. The method of claim 9, wherein a fast parallel algorithm is
used for estimating the updated eigenvalue decomposition.
16. The method of claim 1, further comprising: (e) periodically
generating the covariance matrix based on additional sampled
aperture data; (f) evaluating the eigenvalue decomposition by block
processing each time the covariance matrix is evaluated; (g)
determining a number of jammers within a field of view of the
antenna based on the updated eigenvalue decomposition.
17. A radar processing system, comprising: means for sampling
aperture data received by an antenna array, wherein the sampled
aperture data include data that do not correspond to echo returns
from a beam transmitted by the antenna; means for generating a
covariance matrix using the sampled aperture data; means for
performing an eigenvalue decomposition on the covariance matrix;
means for determining a direction of arrival from which at least
one jammer is transmitting a signal included in the sampled
aperture data, based on the eigenvalue decomposition.
18. The system of claim 17, wherein the direction of arrival is
provided to a tracker.
19. The system of claim 17, wherein the direction determining means
determines an angular location of a jammer that is not within a
mainlobe of the antenna array.
20. The system of claim 17, wherein the direction determining means
determines angular locations of a plurality of jammers that are
located within one or more sidelobes of the antenna array.
21. The system of claim 20, wherein the direction determining means
determines an angular location of an additional jammer that is
within a mainlobe of the antenna array.
22. The system of claim 21, wherein the direction determining means
determines angular locations of a further plurality of jammers that
are within the mainlobe of the antenna array.
23. The system of claim 17, wherein the direction determining means
determines angular locations of a plurality of jammers that are
within a mainlobe of the antenna array.
24. The system of claim 17, wherein the direction determining means
operates while the radar system is operating in surveillance
mode.
25. The system of claim 17, further comprising: means for updating
the covariance matrix repeatedly based on additional sampled
aperture data; means for updating the eigenvalue decomposition each
time the covariance matrix is updated; means for determining a
number of jammers within a field of view of the antenna based on
the updated eigenvalue decomposition.
26. The system of claim 25, wherein the number of jammers within
the field of view is determined based on a number of significant
eigenvalues of the covariance matrix.
27. The system of claim 25, wherein covariance matrix updates are
performed at regular intervals.
28. The system of claim 27, wherein the regular intervals have a
period of between about one millisecond and about ten
milliseconds.
29. The system of claim 25, wherein the direction determining means
uses a super-resolution subspace algorithm to estimate a direction
of arrival of signals from each of the plurality of jammers.
30. The system of claim 25, wherein a sliding window is used to
update the covariance matrix.
31. The system of claim 25, wherein a fast parallel algorithm is
used for estimating the updated eigenvalue decomposition.
32. The system of claim 17, further comprising: (e) means for
periodically generating the covariance matrix based on additional
sampled aperture data; (f) means for evaluating the eigenvalue
decomposition by block processing each time the covariance matrix
is evaluated; (g) means for determining a number of jammers within
a field of view of the antenna based on the updated eigenvalue
decomposition.
Description
FIELD OF THE INVENTION
[0001] The present invention relates to radar systems and methods
generally, and more specifically to methods of detecting a
jammer.
BACKGROUND OF THE INVENTION
[0002] U.S. Pat. No. 5,600,326 to Yu, et al. describes a system for
adaptive beamforming so as to null one mainlobe and multiple
sidelobe jammers. Yu addresses a problem wherein the monopulse
technique for direction of arrival (DOA) estimation failed when
there was sidelobe jamming (SLJ) and/or main lobe jamming (MLJ). If
not effectively countered, electronic jamming prevents successful
radar target detection and tracking. The situation is exacerbated
by introduction of stealth technology to reduce the radar cross
section (RCS) of unfriendly aircraft targets. The frequency
dependence of the RCS encourages use of lower microwave frequency
bands for detection. This leads to large apertures to achieve
angular resolution. Large apertures to achieve small beamwidth
results in interception of more jamming. On the other hand,
constrained apertures lead to wider beamwidth, which implies
interception of more mainlobe jamming.
[0003] Adaptive beamforming techniques have been used to form a
beam having one or more nulls pointing in the direction of one or
more respective jammers. When there is no jamming, Taylor and
Bayliss weightings are typically used for sum beams and difference
beams, respectively, so as to have a narrow mainlobe and low
sidelobes. The quiescent Taylor and Bayliss weightings are designed
for reducing the sidelobes in a practical system. In the presence
of jamming, the weights are adapted so as to form nulls responsive
to the jammers.
[0004] Adaptive receiving arrays for radar maximize the ratio of
antenna gain in a specified scan direction to the total noise in
the output signal. The sum and difference beams at array outputs
are determined by adaptive receiving array techniques, which serve
to null the interference sources. The adaptivity involves using
multipliers to apply an adaptive weight to antenna array signals
furnished at multiplier inputs. Because of the adaptivity, the sum
and difference patterns vary with the external noise field and are
distorted relative to the conventional monopulse sum and difference
beams, which possess even and odd symmetry, respectively, about a
prescribed boresight angle. The adaptive weights for the sum and
difference beams are determined so that the antenna patterns are
distorted. This technique cancels both the mainlobe and sidelobe
jammers but distorts the monopulse ratio.
[0005] Yu et al. describe a sample matrix inverse approach for
jamming cancellation, which effectively forms nulls responsive to
jammers. The covariance matrix is inverted in order to form the
adaptive weighting coefficients. If one of the jammers is within
the mainbeam, a null is formed responsive to the mainlobe jammer
and the mainbeam is distorted. In order to maintain the mainbeam
without distortion, the mainlobe jammer effect is excluded from the
covariance matrix estimate. This may be accomplished by using a
modified covariance matrix in forming the adapted row beamforming
weights, from which information of the mainlobe jammer has been
removed (so there is no null responsive to the mainlobe jammer). Yu
et al. use prefiltering to block the mainlobe jammer.
[0006] Although the matrix inverse approach can generate desired
adaptive weights for pointing nulls toward jammers, this technique
does not output the locations of the jammers. To implement active
countermeasures (such as sending energy at a particular frequency
or band of frequencies to a jammer), it is necessary to know the
location of the jammer. In order to determine the DOA of ajammer
using the prior art techniques, it was necessary to "point" the
receiving antenna array at the jammer, essentially placing the
jammer in the mainlobe. Thus, an improved system is desired.
SUMMARY OF THE INVENTION
[0007] The present invention is a method and system for locating a
radar jammer. Sampled aperture data are received from an antenna
array. The sampled aperture data include data that do not
correspond to echo returns from a beam transmitted by the antenna.
A covariance matrix is generating using the sampled aperture data.
An eigenvalue decomposition is performed on the covariance matrix.
A direction of arrival is determined from which at least one jammer
is transmitting a signal included in the sampled aperture data,
based on the eigenvalue decomposition.
BRIEF DESCRIPTION OF THE DRAWINGS
[0008] FIG. 1 is a diagram showing a plurality of jammers
transmitting signals to a radar antenna array in an exemplary
system according to the present invention.
[0009] FIG. 2 is a block diagram showing exemplary signal
processing performed on the signals received by the antenna array
of FIG. 1.
DETAILED DESCRIPTION
[0010] U.S. Pat. Nos. 5,600,326 to Yu, et al. and 6,087,974 to Yu
are incorporated herein by reference in their entireties, for their
teachings on monopulse radar systems.
[0011] The present invention is a method and system for determining
the angular location of at least one mainlobe or sidelobe jammer.
The preferred embodiment allows determination of the DOA of one or
more mainlobe jammers and/or one or more sidelobe jammers.
Conventional radars locate a jammer by pointing the beam at the
jammer once in a while. In the exemplary embodiment of the present
invention, sampled aperture data, collected during normal target
surveillance and tracking, are used to obtain the jammer DOA and
update the DOA continuously.
[0012] FIG. 1 shows an environment in which the invention may be
practiced. A radar antenna array 100 may receive signals from a
plurality of jammers 120-123. While jammers 120, 121 are located
within the mainlobe 110, jammers 122, 123 are located within
sidelobes. Antenna 100 may be part of a system operating in
surveillance mode. Thus, at different times, the mainlobe may be
pointed in the direction shown by beam 110' or in the direction
shown by beam 110". Note that at different times, jammer 122 may be
in the mainlobe 110' (while jammers 120, 121 and 123 are in
sidelobes) or jammer 123 may be in the mainlobe 110" (while jammers
120-122 are in sidelobes).
[0013] FIG. 2 shows an exemplary system according to the invention.
In the exemplary embodiment, it is not necessary to point the
mainlobe of the antenna 100 in the direction of a jammer to be able
to locate the jammer. The DOA can be determined using sampled
aperture data received by the antenna array 100. The sampled
aperture data include data that do not correspond to echo returns
from a beam transmitted by the antenna 100. Because the jammers
120-123 transmit energy independently of the beam transmitted by
antenna 100, the system can passively collect the sampled aperture
data by collecting data while echo returns from the energy
transmitted by antenna 100 are not being received. Energy from
jammers 120-123, which arrives within the sidelobes, is still
identified. Thus, the sampled aperture data for jammer location
determination can be obtained regardless of whether the radar
system is being operated in surveillance or tracking mode.
[0014] The invention may be practiced using conventional radar
hardware, include an antenna array 100, receivers and
analog-to-digital converters.
[0015] In block 210, a covariance matrix is estimated using the
sampled aperture data. This can be the same covariance matrix
formed to generate the adaptive weighting coefficients. At least
the initial covariance matrix is formed by block processing.
Preferably, updates are estimated to allow rapid, real-time or near
real-time processing, for example, using a sliding window to update
the covariance matrix.
[0016] In block 220, an eigenvalue decomposition is performed on
the covariance matrix. The number of jammers is determined by the
number of significant eigenvalues of the covariance matrix found
while the antenna array 100 is not transmitting. More formal
statistical methods, such as the Akaike Information Criterion (AIC)
or Minimum Description length (MDL) can be used to determine the
number of jammers 120-123 within a field of view of the antenna 100
based on the updated eigenvalue decomposition. Eigenvalue
decomposition may be performed using a conventional block
processing technique. However, in a preferred embodiment, once the
eigenvalue decomposition values are initially determined, a fast
eigenvalue decomposition update (such as an exemplary technique
described below) is performed, to estimate the eigenvalues quickly.
Eigenvalue decomposition is numerically intensive. Thus, a fast
eigenvalue decomposition update technique allows estimation of the
eigenvalues in real-time or near real-time, without interfering
with ongoing surveillance and/or tracking activity.
[0017] Block 230 explicitly determines a direction of arrival from
which at least one jammer 120-123 is transmitting a signal included
in the sampled aperture data, based on the eigenvalue
decomposition. Preferably, a modern super-resolution subspace
algorithm, such as MUSIC, minimum-norm or the like, is used.
[0018] Storage devices 240 and 250 are provided. The storage
devices may be any conventional data storage devices, including
random access memory (RAM), latches, or registers. Each time the
covariance matrix is evaluated (or an update thereof is estimated),
the results are stored in storage device 240, to be used when
updating the covariance matrix on subsequent iterations.
[0019] Similarly, each time the eigenvalue decomposition is
evaluated (or an update thereof is estimated), the results are
stored in storage device 250, to be used when updating the
eigenvalue decomposition on subsequent iterations. Recursive,
cumulative covariance processing allows estimation of the DOA to be
performed in real-time or near real-time with reduced processing
resource consumption.
[0020] In one example, the updates to the covariance matrix are
performed at regular intervals. For example, the regular intervals
may have a period of between about one millisecond and about ten
milliseconds, corresponding to a frequency between about 100 Hz and
1000 Hz. Each time the covariance matrix is updated, the eigenvalue
decomposition is updated, preferably using a fast parallel
decomposition technique. By processing recursively, it is possible
to use the maximum number of samples within the computational
period. Use of sliding window allows the most up-to-date samples to
be used. The covariance processing can be performed continuously in
the background, without interfering with surveillance or tracking
of targets. Further, because the sampled aperture data from any
given direction can be sampled as often as desired,
super-resolution is possible, enabling detection of multiple
jammers within a single range cell.
[0021] Once determined, the direction of arrival may be provided to
a tracker.
[0022] Using the above method, angular locations of one or more
jammers can be determined, including one or more jammers within the
mainlobe of the antenna array 100, and/or one or more jammers
located in sidelobes of the antenna array.
[0023] Because the exemplary method does not require the mainlobe
to be pointed towards a jammer to determine the jammer location, it
is possible to detect and track jammers in normal surveillance
mode.
[0024] In addition to the conventional abilities for target
detection and tracking in the presence of jamming, the exemplary
embodiment makes it possible to add awareness of the DOA of
multiple jammers. This allows the system to monitor and track the
locations of jammers.
[0025] Although the exemplary embodiment described above uses
recursive processing, block processing may alternatively be used.
For block processing, the covariance matrix is periodically
generated based on additional sampled aperture data. The eigenvalue
decomposition is evaluated by block processing each time the
covariance matrix is evaluated. The number of jammers within the
field of view of the antenna is evaluated based on the updated
eigenvalue decomposition. Block processing may be performed every
ten pulses (100 Hz), for example.
FAST EIGENVALUE DECOMPOSITION PROCESSING ALGORITHM
EIGENVALUE DECOMPOSITION
[0026] Eigenspace decompositions are used in solving many signal
processing problems, such as source location estimation,
high-resolution frequency estimation, and beamforming. In each
case, either the eigenvalue decomposition (EVD) of a covariance
matrix or the singular value decomposition (SVD) of a data matrix
is performed. For adaptive applications in a non-stationary
environment, the EVD is updated with the acquisition of new data
and the deletion of old data. This situation arises where the
target sources are moving or the sinusoidal frequencies are varying
with time. For computational efficiency or for real-time
applications, an algorithm is used to update the EVD without
solving the EVD problem from scratch again, i.e., an algorithm that
makes use of the EVD of the original covariance matrix. In
numerical linear algebra, this problem is called the modified
eigenvalue problem. Other examples of modified eigenvalue problems
include extension and restriction problems where the order of the
matrix is modified, corresponding to filter order updating or
downdating.
[0027] One aspect of the exemplary embodiment of the invention is
the use of a fast recursive eigenvalue decomposition method in
processing the radar echo signals in a raid-counting mode.
Formally, these problems are concerned with computing the
eigensystem of Hermitian matrices that have undergone finite
changes that are of restricted rank. The rank values of the
modified part are usually small compared with the rank values of
the original matrices. In these situations, perturbation ideas
relating the eigensystem of the modified and original matrices can
be exploited to derive a computationally efficient algorithm.
[0028] In a time-varying environment, there are two strategies for
updating the covariance matrix. The most common one is the rank-1
update, which involves applying an exponential forgetting window to
the covariance matrix, i.e.,
{circumflex over (R)}=.mu.R+(1-.mu.).alpha..alpha..sup.H (1)
[0029] where R and {circumflex over (R)} correspond to the original
and updated covariance matrices, respectively, .mu. is an
appropriate exponential forgetting factor, and .alpha. is the most
recent data vector. The data vectors correspond to linear
prediction sample vectors in frequency estimation or snapshots in
array processing. The exponential window approach may have several
drawbacks: the influence of old data can last for a very long time,
and the exponential forgetting factor, .mu., must be determined
appropriately.
[0030] Another strategy is to use a sliding window, which is
analogous to short-time signal analysis, where data within the
window is assumed to be stationary. This strategy corresponds to
adding a row to and deleting a row from the data matrix,
simultaneously, or to the following rank-2 covariance matrix update
problem:
{circumflex over (R)}=R+.alpha..alpha..sup.H-.beta..beta..sup.H
(2)
[0031] where .alpha. is the data vector to be included, and .beta.
is the data vector to be deleted. The eigenvalue can be computed
explicitly when the order of the matrix is less than 5 (say, 3 or
4) after deflation. One can also add and delete blocks of data,
leading to the general rank-.kappa. modification problem, and one
can modify equation (2), or the rank-.kappa. modification problem,
by adding weighting factors, thus indicating the relative weight of
the previous covariance matrix and the new data vector in forming
the new covariance matrix estimate.
[0032] Simultaneous data addition and deletion have also been
considered in other signal processing applications such as
recursive least squares problems. Of course, the modification
problem can be solved by applying a rank-1 update as many times as
desired. That approach is computationally involved because
nonlinear searches for eigenvalues have to be carried out for each
rank-1 update, and the procedure is repeated for .kappa. times.
Solving this eigenvalue search problem in an efficient way speeds
up the process. It should also be noted that subtracting data may
lead to ill-conditioning, since the smallest eigenvalue may move
toward zero. Theoretically, this should not happen because only
different segments of data are used for forming the covariance
matrix. Numerically, the ill-conditioning may happen; thus, in
general, data subtraction should be avoided. The choice of updating
schemes depends on specific applications and assumptions in signal
sources. The exemplary embodiment of the present invention provides
an algorithm for recursively updating the EVD of the covariance
matrix when the covariance matrix is under low-rank updating, which
may include data deleting.
[0033] An overview of the modified eigenvalue problem for raid
counting is provided herein and its application to the recursive
updating of the eigen-subspace when the covariance matrix is under
low-rank updating. It includes a theoretical framework in which
rank-.kappa. updates are related to one another. The
spectrum-slicing theorem enables location of the eigenvalue to any
desired accuracy by means of the inertia of a much reduced size
Hermitian matrix whose elements depend on the eigenvalue parameter,
.lambda.. This idea is incorporated and a complete algorithm and
analysis is worked out for updating the eigen-decomposition of the
covariance matrix as it arises in eigenbased techniques for
frequency or angle of arrival estimation and tracking. A novel
algorithm is employed for computing the eigenvector and work out
the deflation procedure for rank-.kappa. updates. The efficient
procedure for computing eigenvalues is a two-step procedure. It
improves the computational complexity from O(N.sup.3) to
O(N.sup.2k). The deflation procedure is important for the
eigen-subspace updating in high-resolution algorithms, as high
filter order leads to multiplicity of noise eigenvalues. If only
the principal eigenvalues and eigenvectors need to be monitored (as
in the PE algorithm), the noise eigenvalue multiplicity
consideration would lead to a highly efficient algorithm. An
analysis of the computational complexity indicates an improvement
of an order of magnitude from O(N.sup.3) to O(N.sup.2k) when
compared with prior known routines. The most efficient EVD routines
involve Householder reduction to tridiagonal form followed by QL
iteration.
[0034] An efficient algorithm for computing the eigensystem of the
modified covariance matrix is described is below. The algorithm
involves a nonlinear search for the eigenvalues that can be done in
parallel. Once the eigenvalues are obtained, the eigenvectors can
be computed explicitly in terms of an intermediate vector, which is
a solution of a deflated homogeneous linear system. The general
procedure, coupled with the multiplicity of noise eigenvalues in
subspace signal processing applications, leads to a highly
efficient algorithm for adaptive estimation or tracking
problems.
[0035] II. Updating the EVD of the Covariance Matrix
[0036] (A) Modified Eigenvalue Problem
[0037] The modified eigenvalue problem is concerned with computing
the eigensystem of the modified Hermitian matrix, given the a
priori knowledge of the eigensystem. of the original matrix. This
specifically concerns the following additive modification:
{circumflex over (R)}=R+E (3)
[0038] where {circumflex over (R)} and R and
R.epsilon.C.sup.N.times.N are the modified and original covariance
matrices and E.epsilon.C.sup.N.times- .N is the additive
modification matrix. This matrix is also Hermitian, and, in
general, is of indefinite nature; i.e., it may have negative
eigenvalues (corresponding to downdating). Assume E is of rank k,
where k is usually much smaller than N. Because E is Hermitian, it
has the following weighted outer product expansion:
E=USU.sup.H (4)
[0039] where U.epsilon.C.sup.N.times.k and
S.epsilon.R.sup.k.times.k is a nonsingular matrix. For example,
equation (4) can be obtained as the eigenvalue and eigenvector
expansion of E, where S is a diagonal matrix with eigenvalues on
the diagonal and U is the corresponding orthonormal eigenvector
matrix. Another example for the decomposition of E, as shown in
equation (4), is expressed directly in terms of the data. In that
case, S has diagonal elements equal to 1 or -1, corresponding to
updating or downdating, respectively, and U is the matrix with the
corresponding data vectors. U is then not orthonormal. A priori
information of the EVD of R is also available as follows:
R=QDQ.sup.H (5)
[0040] where D.epsilon.R.sup.N.times.N, Q.epsilon.C.sup.N.times.N,
D=diag[d.sub.1,d.sub.2 . . . d.sub.N] is the diagonal eigenvalue
matrix, and Q=[q.sub.1, . . . q.sub.N] is the corresponding
eigenvector matrix. Note that Q is an orthonormal matrix, i.e.,
QQ.sup.H=Q.sup.HQ=I (6)
[0041] The problem now is to find the modified eigensystern.
Assuming that .lambda., x are the eigenpairs, of R, the following
expression is obtained:
({circumflex over (R)}-.lambda.I)x=0 (7)
[0042] The eigenvalue can be determined by solving for the zeros
of
det[{circumflex over (R)}-.lambda.I]=0 (8)
[0043] Substituting {circumflex over (R)} from (3) and E from (4)
into (7) gives the following expression:
(R-.lambda.I)x+USU.sup.Hx=0 (9)
[0044] It is convenient to split away the rank-k modification
aspect from the rest of the problem by inducing the following
system of equations, where y=SU.sup.H x and y
.epsilon.=C.sup.k:
(R-.lambda.I)x+Uy=0 (10a)
U.sup.Hx-S.sup.-1y=0 (10b)
[0045] Solving x in terms of y in (10a) and substituting it into
(10b) gives the following:
W(.lambda.)y=0 (11)
where
W(.lambda.)=S.sup.-1+U.sup.H(R-.lambda.I).sup.-1U (12)
[0046] Note that W(.lambda.) can be identified to be the Schur
complement of R-.lambda.I in M(.lambda.), where 1 M ( ) = [ R - I U
U H - S - 1 ] ( 13 )
[0047] W(.lambda.) is also called the Weinstein-Aronszajn (W-A)
matrix, which arises in perturbation problems in Hilbert space
operators. The modified eigenvalues can be obtained from solving
det[W(.lambda.)]=0 rather than from equation (8). In fact, .lambda.
is an eigenvalue of {circumflex over (R)}, i.e.,
.lambda..epsilon..lambda.({circumflex over (R)}) if and only if
.lambda. is a solution of det[W(.lambda.)]=0. This can be derived
easily by using the known Schur equation on equation (13), leading
to 2 det [ R - I ] = det [ M ( ) ] det [ S - 1 ] ( 14 )
[0048] Because S is invertible, det[.left brkt-bot.{circumflex over
(R)}-.lambda.I.right brkt-bot.=0 will imply det[M(.lambda.)]=0.
Also, det[M(.lambda.)] can be expressed as
det[M(.lambda.)]=(-1).sup.kdet[R-.lambda.I]det[W(.lambda.)]
[0049] leading to 3 det [ W ( ) ] = ( - 1 ) k det [ M ( ) ] det [ R
- I ] = ( - 1 ) k i = 1 N ( ^ i - ) i - 1 N ( i - ) ( 15 )
[0050] with {circumflex over
(.lambda.)}.sub.i.epsilon..lambda.({circumfle- x over (R)}) and
.lambda..sub.i.epsilon..lambda.(R). Thus {{circumflex over
(.lambda.)}.sub.i} and {.lambda..sub.i} are the zeros and poles of
the rational polynomial det[W(.lambda.)]. Note that the above
derivation is valid only when the eigenvalues of R and {circumflex
over (R)} are distinct. In fact, R-.lambda.I in equation (12) is
not invertible when .lambda. coincides with one of the eigenvalues
of R. A deflation procedure is prescribed when some of the
eigenvalues of the modified and original matrices are in fact the
same. Note that w(.lambda.)=det[W(.lamb- da.)] is a nonlinear
equation, which is easy to evaluate. W(.lambda.) is a k.times.k
matrix which is of a dimension much smaller than the N.times.N
matrix {circumflex over (R)}-.lambda.I. Moreover, the resolvant
(R-.lambda.I).sup.-1 is easy to compute if the EVD of R is
available, i.e.,
W(.lambda.)=S.sup.-1+U.sup.HQ(D-.lambda.I).sup.-1Q.sup.HU (16)
[0051] where (D-.lambda.I) is a diagonal matrix. The fast algorithm
to be described depends on a spectrum-slicing theorem relating the
eigenvalues of the modified matrix to the eigenvalues of the
original matrix. This theorem enables the search to be localized in
any interval. Moreover, the search can be carried out in parallel,
leading to an efficient implementation.
[0052] Deflation
[0053] Several situations exist for deflating the problem. The
details of two situations are described: (1) U is orthogonal to
q.sub.i, and (2) there are multiplicities of the eigenvalues of the
original covariance matrix. The consideration here is for rank-k
update. Other special situations exist where the problem can be
deflated immediately, such as the existence of zero elements in the
update vector. For the first case, q.sub.i.sup.HU=0; therefore
d.sub.i and q.sub.i remain the eigenpair of the modified matrix
{circumflex over (R)}. For the second case, if .lambda. is an
eigenvalue of multiplicity m with the corresponding set of
eigenvectors Q=[q.sub.i . . . q.sub.m], it is then possible to
perform a Householder transformation on Q, such that the last M-1
eigenvalues are orthogonal to u.sub.1 (where U=[u.sub.1 . . .
u.sub.k]), i.e.,
QH.sub.1=.left brkt-bot.q.sub.1.sup.(1) . . . q.sub.m.sup.(1).right
brkt-bot..DELTA.Q.sup.(1) (17a)
q.sub.i.sup.(1).sup..sup.Hu.sub.1=0 i=2, . . . m (17b)
[0054] Thus, {q.sub.i.sup.(1)}.sub.i=2.sup.m remain the
eigenvectors of {circumflex over (R)} corresponding to eigenvalue
.lambda.. The Householder matrix H.sub.1 is an orthogonal matrix
given by 4 H 1 = 1 - 2 ( b 1 b 1 H ) b 1 H b 1 ( 18 )
[0055] where
b.sub.1=z+.sigma.e.sub.1
[0056] with Q.sup.H u.sub.1=z, .sigma..parallel.z.parallel..sub.2
and e.sub.1=[1,0 . . . 0].sup.H. Now let Q.sub.1=[q.sub.2.sup.(1) .
. . q.sub.m.sup.(1)], and after performing another appropriate
Householder transformation, we obtain
{circumflex over
(Q)}.sup.(1)H.sub.2=[q.sub.2.sup.(2)q.sub.3.sup.(2) . . .
q.sub.m.sup.(2)]q.sub.i.sup.(2)Hu.sub.2=0 i=3, . . . m (19)
[0057] After continuing this procedure k times until
{circumflex over
(Q)}.sup.(k-1)H.sub.k=[q.sub.k.sup.(k)q.sub.k+1.sup.(k). . .
q.sub.m.sup.(k)]q.sub.i.sup.(k)Hu.sub.k=0 i=k+1, . . . m (20)
[0058] then {q.sub.i.sup.(k)}.sub.i=k+1.sup.m are the eigenvectors
of {circumflex over (R)} with .lambda..epsilon..lambda.({circumflex
over (R)}) of multiplicity m-k.
[0059] Spectrum-Slicing Theorem
[0060] Assuming all deflation procedures have been carried out,
i.e., all d.sub.i are distinct and q.sub.i.sup.HU=0, a
computationally efficient algorithm can be used to locate the
eigenvalues to any desired accuracy. The fast algorithm that allows
one to localize the eigenvalue search depends on the following
spectrum-slicing equation:
N.sub.{circumflex over
(R)}(.lambda.)=N.sub.R(.lambda.)+D.sup.+[W(.lambda.- )]-D.sup.+[S]
(21)
[0061] where N.sub.{circumflex over (R)}(.lambda.) and
N.sub.R(.lambda.) are the number of eigenvalues of {circumflex over
(R)} and R less than .lambda., respectively, D.sup.+[W(.lambda..)]
is the positive inertia of W(.lambda.) (i.e., the number of
positive eigenvalues of W(.lambda.) and, likewise, D.sup.+[S] is
the positive inertia of S. The above spectrum-slicing equation (21)
provides information of great computational value. W(.lambda.) is
easy to compute for each value of A. If Q.sup.HU is computed and
stored initially, the dominant cost is on the order of Nk.sup.2
floating-point operations per evaluation. Upon evaluating the W-A
matrix W(.lambda.), one may then compute its inertia efficiently
from an LDL.sup.H or the diagonal pivoting factorization. This
requires k.sup.3 number of operations. The value of
N.sub.R(.lambda.) is available readily, as the EVD of R. D.sup.+[S]
can be determined easily either from the EVD of E (because only a
few eigenvalues are nonzero) or from counting the number of the
data vectors to be added (assuming they do not align with each
other). This needs to be determined only once. The spectrum-slicing
equation provides an easy mechanism for counting the eigenvalues of
{circumflex over (R)} in any given interval. This information
enables one to localize the search in much the same way as Sturm
sequence/bisection techniques are used in polynomial root finding.
The bisection scheme can be carried out until convergence occurs.
The convergence rate is linear. The eigenvalue search algorithm can
be summarized as follows:
[0062] 1. Use the knowledge of the original spectrum and equation
(21) to localize the eigenvalues to disjoint intervals.
[0063] 2. For each iteration step of each eigenvalue search in the
interval (l,u) set 5 = l + u 2 ,
[0064] and test it with N.sub.{circumflex over (R)}(.lambda.)
(equation (21)). Repeat until convergence to desired accuracy.
[0065] Since the distribution of the eigenvalue of the original
matrix is known, the eigenvalues can be localized to disjoint
intervals easily by using equation (21). For the bisection search,
the dominant cost for each iteration is the evaluation of
W(.lambda.) and the LDL.sup.H decomposition for the evaluation of
W(.lambda.).
[0066] This algorithm can be illustrated by considering the
following numerical example. Let the original covariance matrix be
given by R=diag(50, 45, 40, 35, 30, 25, 20, 15, 10, 5). Thus the
eigenvalues are the diagonal elements, and the eigenvectors are the
unit vectors {e.sub.i}.sub.i.sup.10=.sub.1. Now assume the
covariance matrix undergoes the following rank-3 modification given
by
{circumflex over
(R)}=R+u.sub.1u.sub.1.sup.t+u.sub.2u.sub.2.sup.t+u.sub.3u-
.sub.3.sup.t
[0067] where u.sub.1, u.sub.2, and u.sub.3 are randomly generated
data vectors given by [0.3563, -0.2105, -0.3559, -0.3566, 2.1652,
-0.5062, -1.1989, -0.8823, 0.7211, -0.00671.sup.t[-0.5539, -0.4056,
-0.3203, -1.0694, -0.5015, 1.6070, 0.0628, -1.6116, -0.4073,
-0.59501}.sup.t, and [0.6167, -1.1828, 0.3437, -0.3574, -0.4066,
-0.3664, 0.8533, -1.5147, -0.7389, 2.1763].sup.t. Thus
S=diag(1,1,1) and D.sup.+[S]=3. For this example, the resulting
eigenvalues are {circumflex over (.lambda.)}=(51.1439, 47.1839,
40.6324, 36.9239, 36.0696, 27.6072, 22.1148, 20.0423, 11.0808,
8.6086). Now let us illustrate how (21) can be used to locate the
eigenvalues accurately. Start with the interval (20, 25). The
associated W-A matrices are evaluated and the inertias computed.
The counting formulas evaluated at 20+.epsilon. and 25-.epsilon.
(N.sub.R(25-.epsilon.)=4, N.sub.R(20+.epsilon.)=2, .epsilon.=a
small constant depending on accuracy requirement) indicate that
there are two eigenvalues in the interval (20, 25).
[0068] This interval is then split into intervals (20, 22.5) and
(22.5, 25). Evaluating the counting equation indicates that an
eigenvalue has been isolated in each disjoint interval. Bisection
can then be employed to the desired accuracy. Table 1 illustrates
the iteration steps for the interval (20, 21.25) for an accuracy of
10.sup.-3 digits. It converges to 20.0425 in 12 steps. The maximum
number of iterations required for convergence depends on the
interval, (l,u), to be searched and the accuracy requirement,
.epsilon.. This number q can be determined such that
2.sup.q>(u-l)/.epsilon. and is given by Table 2.
1TABLE 1 Bisection Search for Eigenvalues Using the
Spectrum-Slicing Equation Bisection N.sub.R Step Interval Midpoint
N.sub.R(.lambda.) D.sup.+[W(.lambda.)] (.lambda.) 1 (20,21.25)
20.625 4 2 3 2 (20, 20.625) 20.3125 4 2 3 3 (20, 20.3125) 20.1563 4
2 3 4 (20, 20.1563) 20.0782 4 2 3 5 (20, 20.0782) 20.0391 4 1 2 6
(20.0391, 20.0782) 20.0587 4 2 3 7 (20.0391, 20.0587) 20.0489 4 2 3
8 20.0391, 20.0489) 20.044 4 2 3 9 (20.0391, 20.044) 20.0416 4 1 2
10 (20.0416, 20.044) 20.0428 4 2 3 11 (20.0416, 20.0428) 20.0422 4
1 2 12 (20.0422, 20.0428) 20.0425 4 2 3
[0069] When an interval has been found to contain a single
eigenvalue, .lambda., then bisection is a rather slow way to pin
down .lambda., to high accuracy. To speed up the convergence rate,
a preferable strategy would be to switch to a rapidly root-finding
scheme after the roots have been sufficiently localized within
subintervals by using the counting equation. Root-finding schemes,
including variation of secant methods, can be employed. The order
of convergence is 1.618 for secant methods, as against 1 for
bisection. Thus convergence can be accelerated.
[0070] Eigenvector
[0071] Once the eigenvalues are obtained, the eigenvectors can be
evaluated efficiently in two steps. First, the intermediate vector
y can be solved from the k.times.k homogeneous Hermitian system
(11). y can actually be obtained as the byproduct of the LDL.sup.H
decomposition of W(.lambda.); i.e., y is the zero eigenvector of
W(.lambda.) for the convergent eigenvalue .lambda.. The eigenvector
x can then be computed explicitly using equation (10a). This
two-step procedure is much more efficient than solving the original
N.times.N homogeneous system of equation (7) or the conventional
eigenvector solver using the inverse iteration. The explicit
expression x relating y is given by 6 x = - ( R - I ) - 1 U y = - Q
( D - I ) - 1 Q H U y = - i = 1 N q i H U y ( d i - ) q i (22a)
[0072] Normalizing x gives the update eigenvectors {circumflex over
(q)}, i.e., 7 q ^ = x ; x r; 2 (22b)
2 TABLE 2 Number of Number of Subintervals to be Iteration Searched
Steps 10.sup.2 7 10.sup.3 10 10.sup.4 14 10.sup.5 17 10.sup.6
20
[0073] The computational complexity for solving the k.times.k
homogeneous Hermitian system is O(k.sup.3) and the back
substitution of equation (22) involves an order of magnitude
O(N.sup.2k). This is to be contrasted to the Householder reduction
followed by QL Implicit Shift, which requires complexity of order
of magnitude O(N.sup.3). Similarly, the inverse iteration requires
O(N.sup.3). This represents an order of magnitude improvement from
O(N.sup.3) to O(N.sup.2k).
[0074] For completeness, a discussion on the rank-1 and rank-2
updates is provided. Since w(.lambda.) and its derivatives can be
evaluated easily for these two cases, Newton's method can be
employed for the eigenvalue search. Some new results were obtained
as we worked out the complete algorithm, including the deflation
procedure (for high-order rank updates), bounds for the eigenvalues
(derived from the spectrum-slicing theorem instead of DeGroat and
Roberts' generalized interlacing theorem), as well as the explicit
eigenvector computation.
[0075] (B) Rank-1 Modification
[0076] For rank-1 updating modification {circumflex over
(R)}R+.alpha..alpha..sup.H, the modification matrix E can be
identified with E=.alpha..alpha..sup.H, corresponding to the
decomposition equation (4) with U=.alpha. and S=1. W(.lambda.) is a
1.times.1 matrix, which is also equal to the determinant
w(.lambda.), i.e., 8 w ( ) = 1 + H ( R - I ) - 1 = 1 + i = 1 N | q
i H | 2 ( d i - ) ( 23 )
[0077] This is a rational polynomial, with N roots corresponding to
N eigenvalues.
[0078] It is assumed that this is an N.times.N problem for which no
deflation is possible. Thus, all d.sub.i are distinct, and
q.sub.i.sup.Ha.noteq.0. The eigenvalues .lambda..sub.i, of the
rank-1 modified Hermitian matrix satisfy the following interlacing
property:
d.sub.i<.lambda..sub.i<d.sub.i-1i=1,2 . . . N (24a)
[0079] with d.sub.0=d.sub.1+.vertline.a.vertline..sup.2. Therefore,
the search intervals for each .lambda..sub.i can be restricted to
I.sub.i=(d.sub.i, d.sub.i-1). for all i=1,2, . . . N. For
downdating, {circumflex over (R)}=R-.beta..beta..sup.t, the
interlacing equation is given by
d.sub.i+1<.lambda..sub.i<d.sub.i i=1,2, . . . N (24b)
[0080] where d.sub.N+1=d.sub.N-.vertline..beta..vertline..sup.2.
The corresponding search intervals for each .lambda. can be
restricted to I=(d.sub.i+1,d.sub.i). An iterative search technique
can then be used on w(.lambda.) to identify the updated
eigenvalues. The function w(.lambda.) is a monotone increasing
function between its boundary points because 9 w ' ( ) = i = 1 N |
q i H | 2 ( d i - ) 2 > 0 ( 25 )
[0081] Thus, the following Newton's method, safeguarded with
bisection, can be applied in an iterative search to identify the
j-th eigenvalue, .lambda.. For convenience, let us denote the
endpoints of the interval by l=d.sub.j, and u=d.sub.j-1. The full
set of eigenvalues can be solved in parallel by using the following
iterative algorithm for each eigenvalue: 10 * ( k + 1 ) = ( k ) - w
( ( k ) ) w ' ( ( k ) ) ( 26 ) ( k + 1 ) = { * ( k + 1 ) i f * ( k
+ 1 ) I j ( k ) + u 2 i f * ( k + 1 ) u , ( k ) replaces l ( k ) +
l 2 i f ( k + 1 ) l , ( k ) replaces u ( 27 )
[0082] and stopping if
.vertline..lambda..sup.(k+1)-.lambda..sup.(k).vertl-
ine.<.delta..lambda..sup.(k+1), where .delta. is the convergence
threshold. Note that when an iteration goes outside the restricted
zone, since the gradient is greater than zero, this indicates the
direction in which the solution is to be found. Consequently, the
interval can be further restricted as indicated. Newton's method is
guaranteed to converge, and this convergence is quadratic near the
solution. It should be noted that Newton's method is based on a
local linear approximation to the function w(.lambda.). Since
w(.lambda.) is a rational function, it is possible to speed up the
convergence rate by using simple rational polynomials for local
approximations.
[0083] Once the eigenvalues are determined to sufficient accuracy,
the eigenvector can be solved in a procedure described in (24). It
is given explicitly as follows: 11 x = - i = 1 N q i H ( d i - ) q
i (28a) q ^ = x ; x r; 2 (28b)
[0084] Exponential Window Rank-1 Update
[0085] This procedure can be modified to accommodate the
exponential window easily. For the following exponential window,
rank-1 updating modification {circumflex over
(R)}=.mu.R+(1-.mu.)aa.sup.H the Weinstein-Aronszajn matrix
W(.lambda.) is a 1.times.1 matrix, which is also equal to the
determinant w(.lambda.), 12 w ( ) = 1 + ( 1 - ) H ( R - I ) - 1 = 1
+ ( 1 - ) i = 1 N | q i H | 2 ( d i - ) ( 29 )
[0086] The N roots are the updated eigenvalue. They must satisfy
the following modified interlacing property:
.mu.d.sub.i<.lambda..sub.i<.mu.d.sub.i-1 i=1,2, . . . N
(30)
[0087] With ud.sub.0=ud.sub.1+.vertline.a.vertline..sup.2.
Therefore, the search intervals for each .lambda..sub.i can be
restricted to I.sub.i=(.mu.d.sub.i, .mu.d.sub.i-1) for all i=1,2 .
. . N. Assuming eigenvalues do not change dramatically from one
update to another, d.sub.i would be a good initial estimate for
.lambda..sub.i. Once the eigenvalues are determined to sufficient
accuracy, the eigenvector is given explicitly as follows: 13 x ^ i
= - ( 1 - ) k = 1 N ( q k H ) q k ( d k - i ) ( 31 )
[0088] This can be normalized accordingly.
[0089] (C) Rank-2 Modification
[0090] For a rank-2 modification problem, equation (4) can be
expressed in the following: 14 R ^ = R + [ ] [ 1 0 0 - 1 ] [ H H ]
= R + U S U H ( 32 )
[0091] where U=[a.beta.] and S=diag(1, -1). The Weinstein-Aronszajn
matrix W(.lambda.) is now given by
W(.lambda.)=S.sup.-1+U.sup.HQ(D-.lambda.I).sup.-1Q.sup.HU (33)
[0092] Let Q.sup.HU=[yz] and W(.lambda.) can be computed easily
with the following expression for the determinant w(.lambda.) 15 w
( ) = ( 1 + i = 1 N | y i | 2 ( d i - ) ) ( 1 - i = 1 N | z i | 2 (
d i - ) ) + ( i = 1 N y i H z i ( d i - ) ) ( i = 1 N z i H y i ( d
i - ) ) ( 34 )
[0093] Assuming all deflation procedures have been carried out,
i.e., all d.sub.i are distinct and q.sub.i.sup.H
a.noteq.0,q.sub.i.sup.H .beta..noteq.0, simultaneously, the
interlacing property relating the modified eigenvalues to the
original eigenvalues is much more complicated than the rank-1 case.
Combining the interlacing theorems for the rank-1 update equation
(24a) and the rank-1 downdate equation (24b) simultaneously,
provides the following generalized interlacing theorem for the
rank-2 update:
d.sub.i+1<.lambda..sub.i<d.sub.i-1 (35)
[0094] where d.sub.N+1d.sub.N-.vertline..beta..vertline..sup.2 and
d.sub.0=d.sub.1+.vertline.a.vertline..sup.2. That is, the modified
eigenvalues are bounded by the upper and lower eigenvalues, of the
original covariance matrix. In this situation, each interval may
have zero, one, or two eigenvalues. Fortunately, the
spectrum-slicing equation can be used to isolate the eigenvalues in
disjoint intervals. Once the eigenvalues are isolated, Newton's
method can be employed for the eigenvalue search. This nonlinear
search can be implemented in parallel as for the rank-1 case.
[0095] The eigenvector can be determined in two steps. First, an
intermediate vector, y, can be obtained as the solution of the
homogeneous system equation (11), and the eigenvector, x, can be
obtained explicitly in terms of y, as in equation (22). Since k=2,
we can specify the homogeneous system solution y=[1v].sup.H.
Solving W(.lambda.) y=0 gives v 16 v = - ( 1 + a ) b w h e r e a =
1 + i = 1 N | y i | 2 ( d i - ) a n d b = i = 1 N y i H z i ( d i -
) ( 36 )
[0096] Using equation (24) we have the following explicit
expression for the eigenvector 17 x = - i = 1 N y i + v z i ( d i -
) q i a n d q = x ; x r; 2 ( 37 )
[0097] (D) Signal and Noise Subspace Updates with Multiplicity
[0098] In many applications, it is convenient to decompose the
eigenvalues and eigenvectors into signal and noise eigenvalues and
eigenvectors, respectively, i.e., 18 1 2 M M + 1 N 2 ( 38 )
[0099] and
S=[q.sub.1q.sub.2 . . . q.sub.M] (39)
N=[q.sub.M+1 . . . q.sub.N] (40)
[0100] The first M eigenvalues are the signal eigenvalues
corresponding to M complex sinusoids in frequency estimation or M
target sources in array processing. The last N-M eigenvalues
cluster to each other and correspond to the noise level. S and N
are the signal and noise subspace, respectively. For a rank-1
update, the first M+1 eigenvalues and eigenvectors are modified,
and the last N-M-1 eigenpairs stay the same. In order to conform to
the model for M sources, we have the following update equation for
eigenvalues:
{circumflex over (.lambda.)}.sub.1.gtoreq.{circumflex over
(.lambda.)}.sub.2.gtoreq. . . . .gtoreq.{circumflex over
(.lambda.)}.sub.M{circumflex over (.sigma.)}.sup.2 (41)
[0101] where 19 ^ 2 = ^ M + 1 + ( N - M - 1 ) 2 ( N - M ) ( 42
)
[0102] Similarly, for a rank-k update, the first M+k eigenpairs are
modified, and the rest remain the same. We have the following
update equation for the noise eigenvalues according to the model
given by equation (38): 20 ^ 2 = ^ M + 1 + ^ M + 2 + ^ M + k + ( N
- M - k ) 2 ( N - M ) ( 43 )
[0103] If, in fact, there are only M sources, {{circumflex over
(.lambda.)}.sub.M+i}.sub.i=1.sup.k should be close to
.sigma..sup.2. If {{circumflex over
(.lambda.)}.sub.M+1}.sub.i=1.sup.k are not close to .sigma..sup.2
there may be another sinusoid or target. This observation can be
used for detection of new sources.
[0104] (E) Computational Complexity
[0105] Note that rather different techniques are used in the fast
algorithm and in the conventional algorithm for solving the
eigensystem. Basically, the present algorithm involves iterative
search for eigenvalues (using the LDL.sup.H decomposition together
with the spectrum-slicing equation) and then subsequent eigenvector
computation. However, all modern eigensystem solvers involve
similarity transformation for diagonalization. If a priori
knowledge of the EVD of the original covariance matrix were not
available, the similarity transformation technique would be
generally preferred, since it has better computational efficiency
and numerical properties. The present fast algorithm makes use of
the fact that the eigenvalue of the modified matrix is related to
the eigenvalue of the original matrix. Moreover, the eigenvalue can
be isolated to disjoint intervals easily and simultaneously
searched in parallel. It can be located to any desired accuracy by
using a spectrum-slicing equation requiring evaluation of the
inertia of a much-reduced-size (k.times.k) Hermitian matrix. In
this situation, the computational complexity is evaluated and
compared with the complexity of the conventional eigensystem
solver. The results are compared with distinct eigenvalue
situations. Further reduction in computational complexity can be
achieved when the multiplicity of the noise eigenvalues can be
exploited as discussed above in the Signal and Noise Subspace
Updates section.
[0106] The grand strategy of almost all modern eigensystem solvers
is to push the matrix R toward a diagonal form by a sequence of
similarity transformations. If the diagonal form is obtained, then
the eigenvalues are the diagonal elements, and the eigenvectors are
the columns of the product matrices. The most popular algorithms,
such as those used in IMSL, EISPACK, or the Handbook for Automatic
Computation, involve the following two steps: (1) Householder
transformation to reduce the matrix to tridiagonal form and (2) the
QL algorithm with an implicit shift to diagonalize it. In the limit
of large N, the operation count is approximately 21 4 1 3 N 3 .
[0107] For the algorithm described herein, the efficiency depends
on the fact that the eigenvalues can be isolated easily and thus
searched in parallel. For each iteration step, W(.lambda.) is
evaluated and the inertia computed. The overall computational
complexity is summarized as follows:
3 Function No. of Operations Overhead Evaluate S.sup.-1 k Evaluate
Q.sub.1.sup.H = Q.sup.HU N.sup.2k Evaluate {q.sub.liq.sub.1i}.sub.i
=1.sup.N Nk.sup.2 Total N.sup.2k + Nk.sup.2 + k = N.sup.2k
Eigenvalue Evaluate W(.lambda.) Nk.sup.2 LDL.sup.H decomposition
k.sup.3 Total for N.sub.i number of iterations N.sub.i(Nk.sup.2 +
k.sup.3) = N.sub.iNk.sup.2 Eigenvector Backsubstitution N.sup.2(k +
1) - N.sup.2k Total (overhead, eigenvalue, eigenvector) 2N.sup.2k +
N.sub.iNk.sup.2
[0108] Thus, in general, the computational complexity is of order
2N.sup.2k+N.sub.iNk.sup.2. It corresponds to an improvement from 22
4 1 3 N 3 t o 2 N 2 k + N i N k 2 .
[0109] (F) Numerical Properties
[0110] It should be noted that the eigenvalue can be determined to
any desired accuracy by using the spectrum-slicing equation. The
nonlinear search is bisection with a linear convergent rate. The
number of iterations depends on the accuracy requirement, and
details are discussed above in Section II-A. A numerical example is
also included in Section II-A to illustrate how the
spectrum-slicing equation can be used to locate the eigenvalues to
any desired accuracy. The calculation of the updated eigenvectors
depends upon the accurate calculation of the eigenvalues. A
possible drawback of the recursive procedures is the potential
error accumulation from one update to the next. Thus the eigenpairs
should be as accurate as possible to avoid excessive degradation of
the next decomposition, which can be accomplished either from
pairwise Gram-Schmidt or from full-scale Gram-Schmidt
orthogonalization procedures on the derived eigenvectors.
Experimental results indicate that the pairwise Gram-Schmidt
partial orthogonalization at each update seems to control the error
buildup in the recursive rank-1 updating. Another approach is to
refresh the procedure from time to time to avoid any possible
roundoff error buildup.
[0111] III. Adaptive EVD for Frequency or Direction of Arrival
Estimation and Tracking
[0112] This section applies the modified EVD algorithms to the
eigenbased techniques for frequency or angle of arrival estimation
and tracking. Specifically, the adaptive versions of the principal
eigenvector (PE) method of Tufts and Kumersan, the total least
squares (TLS) method of Rahman and Yu, and the MUSIC algorithm are
applied.
[0113] (A) Adaptive PE Method
[0114] In the linear prediction (LP) method for frequency
estimation, the following LP equation is set up: 23 [ x ( 0 ) x ( L
- 1 ) x ( N - L - 1 ) x ( N - 2 ) ] [ c ( L ) c ( 1 ) ] = [ x ( L )
x ( N - 1 ) ] ( 46 )
[0115] or
X.sup.Hc=x (47)
[0116] where L is the prediction order chosen to be
M.ltoreq.L.ltoreq.N-1, M is the number of sinusoids, and X.sup.H is
the data matrix. Note that in the case of array processing for
angle of arrival estimation, each row of the data matrix is a
snapshot of sensor measurement. Premultiplying both sides of
equation (47) by X gives the following covariance matrix normal
equation:
Rc=r (48)
[0117] where R=XX.sup.H and r=Xx. The frequency estimates can be
derived from the linear prediction vector coefficient c. In a
non-stationary environment, it is desirable to update the frequency
estimates by modifying the LP system continuously as more and more
data are available. Specifically, the LP system is modified by
deleting the first k rows and appending another k rows
(corresponding to a rank-2k update) as follows: 24 [ x ( k ) x ( L
+ k - 1 ) x ( N + k - L - 1 ) x ( N + k - 2 ) ] [ c ( L ) c ( 1 ) ]
= [ x ( L + k ) x ( N + k - 1 ) ] ( 49 )
[0118] which leads to the following modification of the covariance
matrix normal equation:
{circumflex over (R)}={circumflex over (r)}(50)
[0119] where 25 R ^ = R + i = 1 k i i H - i = 1 k i i H a n d ( 51
) r ^ = r + i = 1 k i x ( N + i - 1 ) - i = 1 k i x ( L + i - 1 ) (
52 )
[0120] where .alpha..sub.i=[x(N+i-L-1) . . . x(N+i-2)].sup.H and
.beta..sub.i=[x(i-1) . . . x(L+i-2)].sup.H. The PE solution for is
obtained using the following pseudo-rank approximation (assuming
that there are M sources): 26 c ^ = i = 1 M q ^ i H r ^ ^ i q ^ i (
53 )
[0121] where {circumflex over (.lambda.)} and {circumflex over
(q)}.sub.i are the eigenpair solutions at each time instance. The
frequencies can then be obtained by first solving the zeros of the
characteristic polynomials formed by . M of the zeros closest to
the unit circle are then used for determining the sinusoidal
frequency estimates.
[0122] Because only M principal eigenvalues and eigenvectors are
used in equation (53) for solving the LP coefficient vector , it is
desirable to modify the previous algorithm to monitor only the M
principal eigenvalues and eigenvectors instead of the complete set
of N eigenvalues and eigenvectors to facilitate computational
efficiency. In order to implement the update, noise eigenvalues
need to be monitored. It is not necessary to monitor the entire set
of noise eigenvectors, however. Now assume the principal eigenvalue
.lambda..sub.1.gtoreq..lambda..sub.2.gtor- eq. . . .
.gtoreq..lambda..sub.M, the noise eigenvalue .sigma..sup.2 (of
multiplicity N-M), and the signal subspace S=[q.sub.1 . . .
q.sub.M] are available. The rank-2k update algorithm requires the
knowledge of 2k noise eigenvectors, the noise eigenvalues, and the
principal eigenvalues and eigenvectors. The 2k noise eigenvectors
can be obtained in the normalized orthogonal projection of
{.alpha..sub.i} and {.beta..sub.i} on noise subspace N=[q.sub.M+1 .
. . q.sub.N]. These constructions lead to q.sub.j.sup.H
.alpha..sub.i=0 and q.sub.j.sup.H .beta..sub.i=0 for i=1 , . . . k
and j=M+2k, . . . N, and it is not necessary to construct
{q.sub.j}.sub.j=M+2k+1.sup.N. (This procedure corresponds to the
deflation situation and subspace update with multiplicity as
discussed in an earlier section.) .sigma..sup.2 remains as the last
N-M-2k+1 eigenvalue. The algorithm is summarized as follows:
[0123] 1. Construct {q.sub.M+i}.sub.i=1.sup.2k such that they are
orthogonal to {q.sub.i}.sub.i=M+2k+1.sup.N, i.e., 27 q M + j = j -
i = 1 M + j - 1 ( q i H j ) q i j = 1 , k q M + k + j = j - i = 1 M
+ k + j - 1 ( q i H j ) q i j = 1 , k
[0124] 2. Conduct a nonlinear parallel search for {{circumflex over
(.lambda.)}.sub.1, {circumflex over (.lambda.)}.sub.2, . . .
{circumflex over (.lambda.)}.sub.M+2k} using (21). Conduct a
nonlinear parallel search for {{circumflex over (.lambda.)}.sub.1,
{circumflex over (.lambda.)}.sub.2, . . . {circumflex over
(.lambda.)}.sub.M+2k} using (21).
[0125] 3. Update the noise eigenvalue {circumflex over
(.sigma.)}.sup.2 28 ^ 2 = ^ M + 1 + + ^ M + 2 k + ( N - M - 2 k ) 2
( N - M )
[0126] 4. Update the signal subspace eigenvector using (24).
[0127] (B) Adaptive Total Least Squares Method
[0128] The TLS method is a refined and improved method for solving
a linear system of equations when both the data matrix, X.sup.H,
and the observed vector, x, are contaminated by noise. In the LP
method for frequency estimation, both sides of the equation are
contaminated by noise. Thus it is desirable to apply the TLS method
to solve the LP equation. The relative factor for weighting the
noise contribution for both sides of the equation is equal to 1,
because X.sup.H and x are derived from the same noisy samples
{x(i)}. This leads to the following homogeneous system of
equations: 29 [ x ( 0 ) x ( L ) x ( N - L - 1 ) x ( N - 1 ) ] [ c -
1 ] = 0 ( 54 )
[0129] or in terms of the covariance matrix 30 R [ c - 1 ] = 0 ( 55
)
[0130] where R=XX.sup.H and X.sup.H is the data matrix for (54). In
a time-varying environment, we delete k rows and append k rows to
the data matrix X.sup.H in (54) as done in the adaptive PE method
discussed above, leading to the following rank-2k modification of
equation (55): 31 R ^ [ c - 1 ] = 0 ( 56 )
[0131] where 32 R ^ = R + i = 1 k i i H - i = 1 k i H
[0132] and {.alpha..sub.i}, {.beta..sub.i} are the appended and
deleted rows, respectively. The TLS solution is obtained as the
null vector solution of {circumflex over (R)}. In general, for M
sinusoids or M target sources, there are N-M number of minimum
eigenvalues and eigenvectors. Any vector in the noise subspace is a
solution. Out of these infinite number of solutions, one can choose
the following minimum norm solution: 33 c ^ = - k = M + 1 N q ^ k *
L + 1 k = M + 1 N [ q ^ k * ] L + 1 q ^ k ' w h e r e q ^ k = [ q ^
k ' [ q ^ k ] L + 1 ] ( 57 )
[0133] Thus, in the TLS solution, it is the noise subspace
N=[q.sub.M+1 . . . q.sub.N] that must be monitored. Because the
signal eigenvalues and eigenvectors are also used for updating the
noise eigensystem, reduction in computational complexity, as done
in the adaptive PE method, is not achieved.
[0134] (C) Adaptive Music
[0135] In this subsection, the recursive version of a class of
high-resolution algorithms is considered for multiple target angle
estimation or frequency estimation based on the eigenvalue
decomposition of the ensemble-averaged covariance matrix of the
received signal. Consider a system of K moving targets to be
tracked by an array of M sensors. The sensors are linearly
distributed with each sensor separated by a distance d from the
adjacent sensor. For a narrowband signal, the model of the output
of the m-th sensor becomes 34 r m ( t ) = k = 1 K A k ( t ) j2 T k
( t ) ( m - 1 ) d + n m ( t ) m = 1 , 2 , M ( 58 )
[0136] where A.sub.k (t) is the complex amplitude of the k-th
target at time t,T.sub.k(t)=sin {.theta..sub.k(t)} where
.theta..sub.k(t) is the angle of arrival of the k-th target at time
t, and N.sub.m(t) is the m-th sensor noise. Using vector notation,
equation (58) can be written as
r(t)=A(t)s(t)+N(t) (59)
[0137] where 35 r ( t ) = [ r 1 ( t ) r 2 ( t ) r M ( t ) ] s ( t )
= [ s 1 ( t ) s 2 ( t ) s K ( t ) ] n ( t ) = [ n 1 ( t ) n 2 ( t )
n M ( t ) ]
[0138] and the M.times.K direction of arrival (DOA) matrix A(t) is
defined as 36 A ( t ) = [ 1 1 1 j2 T 1 d j2 T 2 d j2 T k d d j2 T 1
( M - 1 ) d j2 T 2 ( M - 1 ) d j2 T k ( M - 1 ) d ]
[0139] The output covariance matrix can then be expressed as
follows:
R(t)=A(t)S(t)A.sup.H(t)+.sigma..sup.2(t)I (60)
[0140] where S(t)=E.left brkt-bot.s(t)s.sup.H(t).right brkt-bot. is
the signal covariance matrix, and .sigma..sup.2(t) is the noise
power. Assuming that K<M, the MUSIC algorithm applied at time t
yields an estimate of the number of targets K, their DOA
{.theta..sub.k(t)}, the signal covariance matrix S(t), and the
noise power .sigma..sup.2(t), by examining the eigenstructure of
the output covariance matrix R(t). R(t) can be estimated from an
ensemble of outer products of snapshots in a sliding window or in
an exponential forgetting window as discussed in Section 1.
[0141] The MUSIC algorithm and its root-finding variations are
briefly reviewed here. Suppose at time t, the estimated covariance
matrix has the following EVD: 37 R = i = 1 N i q i q i H = i = 1 K
i q i q i H + 2 i = k + 1 n q i q i H ( 61 )
[0142] The algorithm depends on the fact that that the noise
subspace E.sub.N=[q.sub.K+1 . . . q.sub.M] is orthogonal to the
signal manifold; i.e., 38 E N H u ( ) = 0 ( 62 )
[0143] where u(.theta.) is the steering vector of the angles to be
searched. The conventional MUSIC algorithm involves searching for
the peaks of the following eigenspectrum: 39 J ( ) = u H ( ) u ( )
u H ( ) E N E N H u ( ) ( 63 )
[0144] To do this, the complete angular interval 40 - 2 2
[0145] is scanned. One can avoid the need for this one-dimensional
scanning by the use of a root-finding approach. This can be
accomplished by, for example, using a known root-MUSIC or minimum
norm algorithm.
[0146] In the root-MUSIC algorithm, 41 j2 T x d
[0147] is replaced by the complex variable z in the eigenspectrum
J(.theta.) defined in (63). Let D(z) denote the resulting
denominator polynomial. The polynomial D(z) can be expressed as the
product of two polynomials, H(z) and H(z.sup.-1), each with real
coefficients. The first polynomial, H(z), has its zero inside or on
the unit circle; K of them will be on (or very close to) the unit
circle and represent the signal zero. The remaining ones represent
extraneous zeros. The zeros of the other polynomial, H(z.sup.-1),
lie on or outside the unit circle, exhibiting inverse symmetry with
respect to the zero of H(z) . The angle estimation is thus
performed by extracting the zeros of the polynomial D(z) and
identifying the signal zeros from the knowledge that they should
lie on the unit circle.
[0148] The minimum norm algorithm is derived by linearly combining
the noise eigenvectors such that:
[0149] 1. The first element of the resulting noise eigenvector is
unity.
[0150] 2. The resulting noise eigenvector lies in the noise
subspace.
[0151] 3. The resulting vector norm is minimum.
[0152] Equation (62) is then modified to 42 A ( ) = 1 H E N E N H u
( ) 1 H E N E N H 1 = 0 ( 64 )
[0153] where .delta.=.sub.1.sup.H=[10 . . . 0]. The angle
estimation problem is then solved by computing the zeros of the
resulting polynomial of equation (64) and identifying the signal
zeros as the K zeros of A(z) that lie on (or very close to) the
unit circle.
[0154] (IV.) Simulation Results
[0155] (A) Numerical Properties
[0156] In this section, the numerical performance of the discussed
algorithms are considered as demonstrated by simulation. The
simulations are performed with a time-varying signal in additive
white noise. Consider the measurement model (Equation 58) for a
scenario where there are three sources (K=3) impinging on a linear
array of 10 sensors (M=10). The signal-to-noise ratio for each
source is 20 dB. The angles are given by
.theta..sub.1(t)=5.degree., .theta..sub.2(t)=22.degree. and
.theta..sub.3(t)=12.degree.. 43 3 ( t ) = 12 .degree. ( K - 1 ) 2
299 .
[0157] In each experiment, 100 updates are carried out. Each EVD
update is obtained by updating the covariance matrix derived from
the data snapshots within a window of length 41.
[0158] As recursive procedures may suffer potential error
accumulation from one update to the next, therefore the sensitivity
of the exemplary algorithm was investigated as a function of the
order of the matrix, and the accuracy requirements for the
eigenvalues searched. Stability and sensitivity tests have been
conducted for this algorithm and comparisons of the performance of
the recursive algorithms with conventional measures.
[0159] The angle estimates for various sizes of matrices (M=7,10)
and eigenvalue search accuracy requirements (tol E-10, E-5) were
compared with the estimates obtained from a conventional routine.
It was observed that the performance is merely a function of the
size of the matrix used, and is not dependent on whether a
conventional eigenvalue decomposition routine or the present
recursive procedure for the eigenvalue decomposition is used, nor
on the accuracy requirements on the eigenvalue search. In fact the
results are practically the same for each case.
[0160] Although the invention has been described in terms of
exemplary embodiments, it is not limited thereto. Rather, the
appended claim should be construed broadly, to include other
variants and embodiments of the invention which may be made by
those skilled in the art without departing from the scope and range
of equivalents of the invention.
* * * * *