U.S. patent number 11,448,774 [Application Number 16/541,730] was granted by the patent office on 2022-09-20 for bayesian geolocation and parameter estimation by retaining channel and state information.
This patent grant is currently assigned to Movano Inc.. The grantee listed for this patent is Movano Inc.. Invention is credited to Matthew C. Bromberg.
United States Patent |
11,448,774 |
Bromberg |
September 20, 2022 |
Bayesian geolocation and parameter estimation by retaining channel
and state information
Abstract
There is provided a method of parameter estimation in a
multi-channel signal environment where a plurality of receive
antennas receive signals from one or more transmitters that
transmit a signal or wave that is reflected from one or more
targets, or receive antennas that receive directly from the
transmitters, whose received signals are processed over multiple
frequencies or channels by a digital receiver. The method comprises
(a) calibrating before the operation of the digital receiver to
determine the free parameters of a mathematical model of a channel
either as the channel model parameters or in the form of tabulated
data; (b) calibrating during the operation of the digital receiver
to determine the channel model; (c) comparing received antennas
array voltages to an analytic or table driven channel model from a
calibrated template without only relying on information lossy
intermediate steps such as delay time of arrival or angle of
arrival measurements; (d) creating a statistical likelihood
function modeling receiver noise to determine model channel
parameters or prior channel uncertainty; (e) saving target
reflector/emitter parameters to be reused for dynamic tracking; and
(f) using Bayesian particle filtering or Maximum Likelihood Methods
to update mixture models for the unknown parameters.
Inventors: |
Bromberg; Matthew C.
(Wenatchee, WA) |
Applicant: |
Name |
City |
State |
Country |
Type |
Movano Inc. |
Pleasanton |
CA |
US |
|
|
Assignee: |
Movano Inc. (Pleasanton,
CA)
|
Family
ID: |
1000006569336 |
Appl.
No.: |
16/541,730 |
Filed: |
August 15, 2019 |
Prior Publication Data
|
|
|
|
Document
Identifier |
Publication Date |
|
US 20200057163 A1 |
Feb 20, 2020 |
|
Related U.S. Patent Documents
|
|
|
|
|
|
|
Application
Number |
Filing Date |
Patent Number |
Issue Date |
|
|
62764814 |
Aug 16, 2018 |
|
|
|
|
Current U.S.
Class: |
1/1 |
Current CPC
Class: |
G01S
19/29 (20130101); H04B 7/0632 (20130101); G06F
17/142 (20130101); G06K 9/6278 (20130101); H04W
4/029 (20180201); G01S 19/256 (20130101); G06V
10/82 (20220101) |
Current International
Class: |
G01S
19/29 (20100101); H04W 4/029 (20180101); G06F
17/14 (20060101); H04B 7/06 (20060101); G06K
9/62 (20220101); G06V 10/82 (20220101); G01S
19/25 (20100101) |
References Cited
[Referenced By]
U.S. Patent Documents
Foreign Patent Documents
|
|
|
|
|
|
|
H 11510981 |
|
Sep 1999 |
|
JP |
|
WO 2004/075577 |
|
Sep 2004 |
|
WO |
|
WO 2013/013169 |
|
Jan 2013 |
|
WO |
|
WO 2016/174679 |
|
Nov 2016 |
|
WO |
|
Other References
Al-Salihi, H. et al., "Bayesian Compressed Sensing-based Channel
Estimation for Massive MIMO Systems," EURASIP Journal on Wireless
Communications and Networking, Dec. 31, 2017, pp. 1-5. cited by
applicant .
PCT International Search Report and Written Opinion, PCT
Application No. PCT/US2019/046737, dated Dec. 5, 2019, 8 pages.
cited by applicant .
PCT International Search Report and Written Opinion, PCT
Application No. PCT/US2019/046741, dated Dec. 5, 2019, 9 pages.
cited by applicant .
Non-Final Office Action, U.S. Appl. No. 16/542,342; (dated Jun. 20,
2022), 17 pgs. cited by applicant.
|
Primary Examiner: Holliday; Jaime M
Attorney, Agent or Firm: Loza & Loza, LLP
Parent Case Text
CROSS REFERENCE TO RELATED INVENTIONS
The present invention claims priority to U.S. Provisional
Application Ser. No. 62/764,814 filed on Aug. 16, 2018, which is
incorporated in its entirety by reference.
Claims
I claim:
1. A method of parameter estimation in a multi-channel signal
environment system where a plurality of receive antennas receive
signals from one or more transmitters that transmit a signal or
wave that is reflected from one or more targets, or receive
antennas that receive directly from the transmitters, whose
received signals are processed over multiple frequencies or
channels by a digital receiver and one or more processors whose
processing comprises the steps of: calibrating before the operation
of the digital receiver to determine free parameters of a
mathematical model of a channel either as channel model parameters
or in the form of tabulated data; calibrating during the operation
of the digital receiver to determine a channel model; comparing
received antennas array voltages to an analytic or table driven
channel model from a calibrated template without only relying on
information lossy intermediate steps such as delay time of arrival
or angle of arrival measurements; creating a statistical likelihood
function modeling receiver noise to determine model channel
parameters or prior channel uncertainty; and saving target
reflector/emitter parameters to be reused for dynamic tracking; and
using Bayesian particle filtering or Maximum Likelihood Methods to
update mixture models for unknown parameters.
2. The method of parameter estimation of claim 1, further including
using Bayesian detection or other Statistical Signal Processing
Techniques for the estimation of channel parameters such as
location parameters, shape parameters and reflector electromagnetic
parameters.
3. The method of parameter estimation of claim 1, wherein the
transmitted signal is a frequency stepped radar.
4. The method of parameter estimation of claim 1, wherein a
statistical likelihood function is used to determine target type or
target position.
5. The method of parameter estimation of claim 1, wherein a channel
template for each target type is used to further classify the
target type or estimate a target position.
6. The method of parameter estimation of claim 1, wherein static
direct path clutter is removed via a cancellation algorithm.
7. The method of parameter estimation of claim 1, that blindly
calculates an unknown gain and phase constant .beta.q over each
channel observation q, to absorb bulk phase and gain changes
between observations.
8. The method of parameter estimation of claim 1, wherein the
statistical likelihood function takes the form of a Cauchy quotient
and is therefore a member of the numerical range of a parameterized
matrix as exhibited in Equation (14).
9. The method of parameter estimation of claim 1, wherein a target
classification is determined from the likelihood function using the
posterior probability of detection via Bayes theorem.
10. The method of parameter estimation of claim 9, wherein the
target classification is determined from the statistical likelihood
function after blindly estimating phase terms such as time delays
of arrival.
11. The method of parameter estimation of claim 1, wherein a target
classification is determined from the statistical likelihood
function after blindly estimating any position dependent parameters
during the classification process itself.
12. The method of parameter estimation of claim 11, wherein the
target classification is determined from the statistical likelihood
function after blindly estimating phase delays using fast Fourier
transform (FFT) processing.
13. The method of parameter estimation of claim 11, wherein the
target classification or position parameters are determined from
the statistical likelihood function that assumed scale parameters
.beta.q=.beta., are independent of a collection index and can
therefore be solved using the generalized eigenvalue problem in
Equation (8).
14. The method of parameter estimation of claim 1, wherein DSP
processing hardware is configured to enable parallel statistical
likelihood function calculations wherein the same instruction
across all processing resources.
15. The method of parameter estimation of claim 1, wherein DSP
hardware has built in support for parallel reduction for primitive
associative operations like addition, maximization, minimization or
multiplication to enable the parallel computation of the
statistical likelihood function.
16. The method of parameter estimation of claim 1, wherein multiple
transmit and receive RF chains are packed into an radio frequency
(RF) integrated circuit (IC) with clocks controlled by a phased
locked loop.
17. The method of parameter estimation of claim 1, wherein the
signals from each antenna are down converted and digitized by a
bank of analog to digital converters, where they are further
decimated to obtain a sample of the channel for each antenna and
frequency in use.
18. The method of parameter estimation of claim 17, further
includes blindly calculating an unknown gain and phase constant
.beta.q over each channel observation q, to absorb bulk phase and
gain changes between observations.
19. The method of parameter estimation of claim 17, wherein DSP
processing hardware is configured to enable parallel statistical
likelihood function calculations, yet can share the same
instruction across all processing resources.
20. The method of parameter estimation of claim 17, wherein DSP
hardware has built in support for parallel reduction for primitive
associative operations like addition, maximization, minimization or
multiplication to enable the parallel computation of the likelihood
function.
21. The method of parameter estimation of claim 17, wherein
multiple transmit and receive radio frequency (RF) chains are
packed into an RF integrated circuit (IC), with clocks controlled
by a phased locked loop.
22. The method of parameter estimation of claim 17, wherein the
signals from each antenna are down-converted and digitized by a
bank of analog to digital converters wherein the signals are
further decimated to obtain a sample of the channel for each
antenna and frequency in use.
23. The method of parameter estimation of claim 17, wherein the
statistical likelihood function is itself used to further refine
presumed positions of the transmitter and receiver to improve a
calibration channel model.
24. The method of parameter estimation of claim 1, further
interpolating a wave function using basis functions that are
reflected through objects of predetermined shape or presumed
electrical properties.
25. The method of parameter estimation of claim 24, wherein
interpolating uses a likelihood function to compute a fixed set of
channel vectors .alpha., modulating the basis functions, that can
be used to interpolate the channel at an arbitrary target
position.
26. The method of parameter estimation of claim 24, further
including computing separate channel vectors .alpha..sub.c for each
target type c.
27. The method of parameter estimation of claim 24, further
including jointly computing the phase constants .beta.q and the
channel vectors .alpha..sub.c by alternating directions
optimization.
28. The method of parameter estimation of claim 24, further
includes assuming an inner product of the basis functions is
independent of position resulting in an eigenvalue problem (10)
that is used to learn the channel vectors during calibration.
29. The method of parameter estimation of claim 24, further
including assuming scale parameters .beta.q=.beta., are independent
of a collection index, and can therefore be solved using a
generalized eigenvalue problem in Equation (8).
30. The method of parameter estimation of claim 24, further
calibrating a blind classifier by blindly estimating any position
dependent parameters during a classification process, without
necessarily using known target or array transceiver positions.
31. The method of parameter estimation of claim 24, further
calibrating a blind classifier by after blindly estimating phase
terms such as time delays during the classification process itself
without necessarily using known target or array transceiver
positions.
32. The method of parameter estimation of claim 24, further
calibrating a blind classifier by blindly estimating phase delays
using fast Fourier transform (FFT) processing, during the
classification process itself, without necessarily using known
target or array transceiver positions.
Description
BACKGROUND OF THE INVENTION
The problem of parameter estimation for electromagnetic signals
remains a central one for geolocation, tracking, pattern
recognition, medical imaging, threat detection and for the
detection of gestures and other time varying behavior. Most
geolocation techniques used rely on phase and angle measurements or
line of bearing measurements which are known to be sub-optimal.
Recently, however, techniques have arisen which do not take the
information lossy step of only saving time of arrival or angle of
arrival, but instead consider the entire channel. These techniques
have been extended for the pure geolocation application in the
presence of interference. However, there are additional missing
considerations for these techniques to achieve their full
potential.
First and foremost an adequate model of the channels in question as
well as an understanding of the model errors and non-idealities are
needed. Indeed previous techniques are potentially very sensitive
to calibration errors, and furthermore it does not address near
field channels which need to be modeled in order to handle near
field applications, such as gesture modeling, shape modeling,
medical imaging and security applications. These simple channel
models are also inadequate to handle radar, especially in the near
field.
There is also a need to marry the extensive literature and novel
techniques in the field of Bayesian estimation, since this allows
one to naturally obtain error estimates and it is well suited for
applying machine learning models for pattern recognition. In
particular one or more embodiments of the invention look at models
of Gaussian mixtures and solve them using techniques similar to
Gaussian Sum particle filters. Gaussian sum filters have been used
in part for the tracking problem, but using GPS to obtain
geolocation. Techniques provided lieleni will apply to both the
near and far field and be more accurate due to the removal of
information lossy steps.
Adding some novel techniques for obtaining solutions to
electromagnetic boundary problems, along with some cutting edge
techniques in Bayesian estimation, and by retaining channel state
information directly, one or more embodiments are directed to
obtain a novel approach to the estimation and detection problem for
multiple antennas/transmitters.
SUMMARY OF THE INVENTION
One or more emboidments of the invention are directed to the use of
multi-antenna radar to identify and classify objects of interest as
depicted in FIG. 1. However, the technology can apply to any device
that transmits and receives over multiple antennas or that receives
over diverse positions in space or transmits over diverse
transmitters. An example of the latter would be a GPS application
as shown in FIG. 2.
The baseline transceiver block diagram is shown in FIG. 3. The
transceiver consists of a clock and tuner that phase synchronizes
both the transmitter and the receiver. The transmit waveform is
typically a stepped frequency tuner, that transmits pure tones over
multiple frequencies, sweeping over those frequencies at a fixed
rate. The receiver receives those tones after the signal is
possibly reflected from an external object. For other applications
the transmitter is not present and the receiver is simply
channelized.
The direct path transmit waveform is isolated from the receivers
via RF shielding, and delay filtering if needed (described below).
If a change is detected in the environment, a position estimate is
formed using Bayesian estimation with a likelihood function formed
using a channel model that models the electromagnetic scattering
over the presumed target type. The estimation is essentially
performed using hypothesis testing on the target position, coupled
with prior information about previous targets.
Both the target type (shape or dielectric constants etc.) are
determined by computing the Bayesian posterior probability of that
type category. The Bayesian estimation is described below in
Parameter Estimation and Tracking. The channel model uses data
obtained through careful calibration of both the antenna arrays,
the antenna noise models, and the scattering type. The channel
model is described in below in System Model and the calibration
procedure is described in Calibration. A more detailed step by step
description of the calibration process is illustrated in FIG. 4.
Each target type must go through an independent calibration
procedure, since the target type will determine the basis function
coefficients.
Numerous other advantages and features of the invention will become
readily apparent from the following detailed description of the
invention and the embodiments thereof, from the claims, and from
the accompanying drawings.
BRIEF DESCRIPTION OF THE FIGURES
The patent or application file contains at least one drawing
executed in color. Copies of this patent or patent application
publication with color drawing(s) will be provided by the Office
upon request and payment of the necessary fee. A fuller
understanding of the foregoing may be had by reference to the
accompanying drawings, wherein:
FIG. 1 is a radar transceiver;
FIG. 2 is a GPS receiver from multiple satellite transmitters;
FIG. 3 is a transceiver block diagram;
FIG. 4 is a calibration procedure;
FIG. 5 is a delay impulse;
FIG. 6 is a Time delay filtering;
FIG. 7 is a hand wave gesture detection algorithm;
FIG. 8 is a multi-layer neural network;
FIG. 9 is a learning network based on higher order SVDs
(HOSVD);
FIG. 10 is a integration box for boundary value calculations;
FIG. 11 shows the HOSVD singular values of a Green's function;
FIG. 12 is a locus plot of maximal singular functions; and
FIG. 13 shows a dipole image through a reflector.
DESCRIPTION OF THE INVENTION
While the invention is susceptible to embodiments in many different
forms, there are shown in the drawings and will be described in
detail herein the preferred embodiments of the present invention.
It should be understood, however, that the present disclosure is to
be considered an exemplification of the principles of the invention
and is not intended to limit the spirit or scope of the invention
and/or claims of the embodiments illustrated.
System Model
The invention proposes in one or move embodiments a general
receiver model that can be utilized for localization, object
classification or even medical imaging. The model focuses on a set
of views of the environment that depend on a parameter vector of
interest p. Consider a diverse set of K observations of an emitter
(or passive scatterer) on one or more multi-antenna
receivers/sensors. The observations are indexed by k for k=1:K. The
"receivers" may in fact be simply a view of the signal of interest
at K different frequencies, or it may be K receivers at different
spatial locations, possibly separated by time, it may be a single
receiver, receiving a signal from K transmitters or it may be a
view of the signal after correlating over K different basis
functions or code words.
The received signal by the k'th receiver (or during the k'th
observation) is modeled by,
.alpha. ##EQU00001## where x.sub.k(t) is the M.sub.k.times.1
received waveform vector at sensor k and time t, i.sub.k(t) is an
additive interference plus noise process for sensor k, .sub.kq
(.alpha..sub.kq,p,s.sub.kq(t)) is a vector channel process that
modulates the transmitted signal s.sub.kq(t). Sensor k receives
signals from possibly Q different emitters or from Q different
reflectors. The channel process depends on a set of possibly random
unknown parameters unique to sensor k and emitter Q, .alpha..sub.kq
and a global parameter vector p.sub.q, which is of primary interest
to our receiver processing algorithms. The p.sub.q parameter might
consist of geolocation coordinates, target velocity, X-ray mass
absorption coefficients, or other signal invariants of
interest.
If it is assumed that each view of the environment is statistically
independent (though this requirement can be relaxed), one obtains
the following likelihood function,
.function..intg..alpha.
.times..alpha..times..alpha..times..times..alpha. ##EQU00002##
where P(x|p) is the probability of observing the total view column
vector x=[x.sub.1.sup.T,x.sub.2.sup.T, . . . ,x.sub.K.sup.T].sup.T,
given p=[p.sub.1.sup.Tp.sub.2.sup.T, . . . p.sub.Q.sup.T].sup.T,
f.sub.k (x) is the probability density function for the
interference process i.sub.k, and g.sub.k (.alpha..sub.kq) is the
prior distribution of the noisy, free parameters of the channel
process .sub.k ( ).
Linear Channel Model
Often of great interest is the narrow band channel model, which
treats the channel process as a simple complex multiply. This
approximation is valid for linear time invariant systems, whose
signal of interest is either narrow band, or has been channelized
into several narrow band channels. In this case, in part the k
index can serve as a channel index. Equation (1) is written as,
.alpha..times. ##EQU00003## The notation can be further simplified
by writing this in the time sampled digital domain as,
x.sub.k(n)=H.sub.ks.sub.k(n)+i.sub.k(n),
H.sub.k.ident.[.sub.k1(.alpha..sub.k1,p.sub.1),.sub.k2(.alpha..sub.k2,p.s-
ub.2) . . . ,.sub.kQ(.alpha..sub.kQ,p.sub.Q)] where n is the time
sample index, H.sub.k is the M.sub.k.times.Q channel matrix as a
function of the nuisance parameter vector
.alpha..sub.k.ident.[.alpha..sub.k1.sup.T, .alpha..sub.k2.sup.T . .
. .alpha..sub.kQ.sup.T] and parameter of interest vector
p.ident.[p.sub.1.sup.T,p.sub.2.sup.T, . . . p.sub.Q.sup.T].sup.T,
for all Q emitters/reflectors. Also the emitter signal waveforms
are packed into the Q.times.1 signal vector
s.sub.k(n).ident.[s.sub.k1(n),s.sub.k2(n) . . .
,s.sub.kQ(n)].sup.T.
Colored Gaussian Interference
Some processing simplifications can be made by subsuming all
unmodeled interfering waveforms into the single interference vector
i.sub.k (n) and further presuming that this interference is complex
Gaussian, yet not white. In particular it can be presumed the
interference takes the form of a zero mean multivariate, complex
circularly symmetric Gaussian vector,
p.sub.i(i.sub.k|R.sub.ii)=|.pi.R.sub.i.sub.k.sub.i.sub.k|.sup.-1
exp(-i.sup.HR.sub.i.sub.k.sub.i.sub.k.sup.-1i). The zero mean
Gaussian normal distribution is referred to by, N (0,
R.sub.i.sub.k.sub.i.sub.k), and one can also define the complex
normal distribution function as, N(i,R).ident.|.pi.R|.sup.-1
exp(-i.sup.HR.sup.-1i). This allows the likelihood function to be
written as,
.function..pi..times..times..times..function..function..times..times..fun-
ction..times. ##EQU00004## where Tr( ) is the matrix trace operator
and o(a).ident.aa.sup.H is the outer product operator.
We now perform some manipulations to expose the sufficient
statistics for the probability function as follows:
.function..pi..times..times..times..function..times..times..times..functi-
on..times..pi..times..times..times..times..times..times..times..times..tim-
es..times..times..times. ##EQU00005## where the time averaged cross
and auto correlation matrices are defined by
.ident..times..function..times. ##EQU00006## and where etr
(A).ident.exp (Tr(A)). A matrix version of completing the square
provides a way to identify sufficient statistics for this problem.
We write:
p(x.sub.k|H.sub.k,s.sub.k)=|.pi.R.sub.i.sub.k.sub.i.sub.k|.sup.-N
etr(-NR.sub.i.sub.k.sub.i.sub.k.sup.-1({circumflex over
(R)}.sub.x.sub.k.sub.x.sub.k-{circumflex over
(R)}.sub.x.sub.k.sub.s.sub.k{circumflex over
(R)}.sub.s.sub.k.sub.s.sub.k.sup.-1{circumflex over
(R)}.sub.x.sub.k.sub.s.sub.k.sup.H+ . . . (H.sub.k-{circumflex over
(R)}.sub.x.sub.k.sub.s.sub.k{circumflex over
(R)}.sub.s.sub.k.sub.s.sub.k.sup.-1){circumflex over
(R)}.sub.s.sub.k.sub.s.sub.k(H.sub.k-{circumflex over
(R)}.sub.x.sub.k.sub.s.sub.k{circumflex over
(R)}.sub.s.sub.k.sub.s.sub.k.sup.-1).sup.H)). Note that
H.sub.k.ident.{circumflex over
(R)}.sub.x.sub.k.sub.s.sub.k{circumflex over
(R)}.sub.s.sub.k.sub.s.sub.k.sup.-1 is a least squares
approximation of the channel matrix, and that {circumflex over
(R)}.sub.i.sub.k.sub.i.sub.k.ident.{circumflex over
(R)}.sub.x.sub.k.sub.x.sub.k-{circumflex over
(R)}.sub.x.sub.k.sub.s.sub.k{circumflex over
(R)}.sub.s.sub.k.sub.s.sub.k.sup.-1{circumflex over
(R)}.sub.x.sub.k.sub.s.sub.k.sup.H is an estimate of the
interference covariance matrix. This allows one to write,
p(x.sub.k|H.sub.k,s.sub.k)=.eta.etr(-NR.sub.i.sub.k.sub.i.sub.k.sup.-1((H-
.sub.k-H.sub.k){circumflex over
(R)}.sub.s.sub.k.sub.s.sub.k(H.sub.k-H.sub.k).sup.H))
.eta..ident.|.pi.R.sub.i.sub.k.sub.i.sub.k|.sup.-N
etr(-NR.sub.i.sub.k.sub.i.sub.k.sup.-1{circumflex over
(R)}.sub.i.sub.k.sub.i.sub.k). (4) This makes
H.sub.k.ident.{circumflex over
(R)}.sub.x.sub.s.sub.s.sub.k{circumflex over
(R)}.sub.s.sub.k.sub.s.sub.k.sup.-1 a sufficient statistic given
R.sub.i.sub.k.sub.i.sub.k and s.sub.k.
We can write this in the whitened form,
p(x.sub.k|H.sub.k,s.sub.k)=.eta.etr(-N(H.sub.sk-{circumflex over
({tilde over (H)})}.sub.sk).sup.H({tilde over
(H)}.sub.sk-{circumflex over ({tilde over (H)})}.sub.sk)) (5) where
{tilde over
(H)}.sub.sk.ident.R.sub.i.sub.k.sup.-HH.sub.k{circumflex over
(R)}.sub.s.sub.k.sup.H is the interference whitened and signal
power normalized channel matrix for view k, and where {dot over
({tilde over (H)})}.sub.sk.ident.R.sub.i.sub.k.sup.-H{circumflex
over (R)}.sub.x.sub.k.sub.s.sub.k{circumflex over
(R)}.sub.s.sub.k.sup.-1, is an estimate of the interference
whitened and signal power normalized channel matrix. The matrix
R.sub.i.sub.k is the upper triangular Cholesky factor of
R.sub.i.sub.k.sub.i.sub.k with R.sub.i.sub.k.sup.-H being the
Hermitian transpose of the inverse of that matrix. Similarly
{circumflex over (R)}.sub.s.sub.k.sup.H is the Hermitian transpose
of the Cholesky factor of {circumflex over
(R)}.sub.s.sub.k.sub.s.sub.k.
The interference matrix and it's Cholesky factors can be estimated
independently, either by observing the system in a quasi-stationary
state, or by canceling out all known signals of interest.
This result is dependent on s.sub.kq(n) being a sufficiently wide
band signal to make {circumflex over (R)}.sub.s.sub.k.sub.s.sub.k a
full rank matrix. We must also cover the very important case of
{circumflex over (R)}.sub.s.sub.k.sub.s.sub.k being a rank one
matrix, since an important special case is where we have either a
channelized system, with channel index k, or a stepped-frequency
radar. For this case, it is easiest to presume that H.sub.k is
actually a rank-1 vector modeled as the sum of it's multipath
components,
.ident..times..alpha. ##EQU00007## Here we assume that the position
parameters p is shared by all reflected multipath components.
Narrow Band Channels
The important case wherein s.sub.kq(n) is a narrow band signal,
such as a sinusoid requires special treatment. This case would
arise, for example when the signal is a multi-tone signal common in
stepped frequency radar, or a known signal type that has been
channelized.
If one specializes to the case wherein k represents a frequency
f.sub.k and q is associated with a path delay .tau..sub.q, we can
approximate the signal of interest as,
s.sub.kq(n).apprxeq.s.sub.k(n)e.sup.-2.pi.jf.sup.k.sup..tau..sub.q,
where s.sub.k (n) is a narrow band signal centered at f.sub.k.
{circumflex over (R)}.sub.x.sub.k.sub.s.sub.k={circumflex over
(r)}.sub.x.sub.k.sub.s.sub.kv.sub.k.sup.H, where
.times..times..times..times..ident..times..pi..times..times..tau..times..-
pi..times..times..tau..times..times..times..pi..times..times..tau.
##EQU00008## Similarly we have, {circumflex over
(R)}.sub.s.sub.k.sub.s.sub.k.ident.{circumflex over
(r)}.sub.s.sub.k.sub.s.sub.kv.sub.kv.sub.k.sup.H.
This allows Equation (3) to be written as,
.function..pi..times..times..times..times..times..times..times..times..ti-
mes..times..times..times..times..times..times..times..pi..times..times..ti-
mes..times..times..times..times..times..times..ident..times..times..times.-
.times..times..ident..times..times..times..times..times.
##EQU00009## Similar to Equation (5), one can write the whitened
version of the likelihood function for the narrow band channel as,
p(x.sub.k|H.sub.k,s.sub.k)=.eta. exp(-N{circumflex over
(r)}.sub.s.sub.k.sub.s.sub.k.parallel.{tilde over (h)}.sub.k-{tilde
over (h)}.sub.k.parallel..sup.2) (7)
.eta..ident.|.pi.R.sub.i.sub.k.sub.i.sub.k|.sup.-N
etr(-NR.sub.i.sub.k.sub.i.sub.k.sup.-1{circumflex over
(R)}.sub.i.sub.k.sub.i.sub.k), where {tilde over
(h)}.sub.k.ident.R.sub.i.sub.k.sup.-Hh.sub.k and {tilde over
(h)}.sub.k.ident.R.sub.i.sub.k.sup.-Hh.sub.k.
Aperture Noise Model
For proper estimation of channel parameters it is critical to have
an accurate noise model of the aperture vector .sub.kq
(.alpha..sub.kq,p.sub.q). This is because some forms of
intereference whitening can cause the true aperture vector to be
canceled in the estimator derived from a naive implementation of
maximum likelihood estimation.
For simplicity we refer to the parameterized aperture vector as
h=.sub.kq (.alpha..sub.kq, p.sub.q). If we subsume all the other
reflectors/receivers/transmitters associated with the other q
indices into the interference vector, from Equation (2) into an
interference vector, a narrow band channel receiver model can be
written for the aperture vector as, x(n)=hs(n)+i(n). From this we
can consider several different aperture noise models.
Additive Gaussian Noise
This model simply adds Gaussian noise to the true aperture h.sub.t,
h=h.sub.t+e, where e is colored Gaussian complex noise. In general
e may be dependent on the unknown parameters .alpha..sub.kq,
p.sub.q. It can be useful to consider a piecewise approximation of
the e noise parameters for use in Bayesian estimation. In
particular we could presume that the covariance matrices are
dependent on an index that is unique to a given region or
assumption about the nuisance parameters.
The whitened least squares estimators suggested by Equation (5) and
Equation (7) are conducive to analyzing aperture noise models under
the assumption that they are Gaussian.
Multiplicative Aperture Noise
The channel may be subject to phase noise due to oscillator
imperfections. This is often modeled as multiplicative noise. In
this case we would model the channel as,
h=h.sub.t.circle-w/dot.(1+e), where .circle-w/dot. represents
elementwise multiplication and e is colored Gaussian complex
noise.
One approximation subsumes this noise model into the simpler noise
model of the Additive Gaussian Noise, with e replaced with
h.sub.t.circle-w/dot.e, another treats this model by an aperture
square gain term added to the diagonal of the interference
covariance matrix, further multiplied by relative phase noise
power.
A more constrained approach to aperture noise allows a noise
multiplier for each spread frequency. So for frequency k we model
the receiver channel as an M.sub.k.times.1 vector, {tilde over
(h)}.sub.k.alpha..sub.k, where {tilde over (h)}.sub.k is a possibly
whitened component of the channel seen at frequency k across all
receivers. The receiver model can be written as, {tilde over
(x)}.sub.k={tilde over (h)}.sub.k.alpha..sub.k+e.sub.k
.alpha..sub.k=1+.epsilon..sub.k. Here, {tilde over (x)}.sub.k is
the whitened received data, in the form of an estimated channel
response, playing the role of {tilde over (h)}.sub.k in Equation
(7).
Assuming a complex Gaussian noise model we have,
.function..alpha..pi..times..sigma..times..times..beta..times..function..-
sigma..times..times..alpha..beta..times..alpha. ##EQU00010## Here
.sigma..sup.2 is the variance of the additive white Gaussian noise
(post whitening) and .beta..sub.k.sup.2 is the variance that
controls the deviation of the per-frequency channel gain from
unity. By varying .beta..sub.k.sup.2 we can express our confidence
in the fidelity of the channel model. A large .beta..sub.k.sup.2
suggests that there is a lot of frequency dependent model error in
the channel, perhaps due to phase noise arising from a stepped
frequency tuner.
We can now integrate away the nuisance parameter representing the
gain .alpha..sub.k to obtain, P({tilde over (x)}.sub.k|{tilde over
(h)}.sub.k)=det(.pi.R.sub.x.sub.k.sub.x.sub.k).sup.-1 exp(-({tilde
over (x)}.sub.k-{tilde over
(h)}.sub.k).sup.HR.sub.x.sub.k.sub.x.sub.k.sup.-1({tilde over
(x)}.sub.k-{tilde over (h)}.sub.k)) (9) where
R.sub.x.sub.k.sub.x.sub.k.ident..sigma..sup.2I+{tilde over
(h)}.sub.k{tilde over (h)}.sub.k.sup.H.beta..sub.k.sup.2. The
posterior expectation for .alpha..sub.k and it's maximum likelihood
estimate is given by,
.alpha..beta..times..times..sigma..beta..times..times..sigma.
##EQU00011## It can be shown that to evaluate the argument in the
exponential of Equation (9), it suffices to substitute {circumflex
over (.alpha.)}.sub.k into .alpha..sub.k in Equation (8). Also the
determinant is easily computed as,
det(.pi.R.sub.x.sub.k.sub.x.sub.k)=.pi..sup.M.sup.k(.sigma..sup.2+.beta..-
sub.k.sup.2.parallel.{tilde over (h)}.sub.k.parallel..sup.2).
In the case that the multiplicative noise is applied over both
antenna number and frequency number, we have the receiver channel
vector modeled as, {tilde over (x)}.sub.mk={tilde over
(h)}.sub.mk.alpha..sub.mk+e.sub.mk, where m is the antenna index
and k is the frequency index. This model leads to the conditional
distribution of, P({tilde over (x)}|{tilde over
(h)}.sub.k)=det(.pi.R.sub.x.sub.k.sub.x.sub.k).sup.-1 exp(-({tilde
over (x)}.sub.k-{tilde over
(h)}.sub.k).sup.HR.sub.x.sub.k.sub.x.sub.k.sup.1({tilde over
(x)}.sub.k-{tilde over (h)}.sub.k)), where
R.sub.x.sub.k.sub.x.sub.k.ident..sigma..sup.2I+(|{tilde over
(h)}.sub.k|.sup.2)(.beta..sub.k.sup.2), where (|{tilde over
(h)}.sub.k|.sup.2) is an M.times.M diagonal matrix, whose m'th
diagonal element is |{tilde over (h)}.sub.mk|.sup.2 and
(.beta..sub.k.sup.2) is a diagonal matrix whose m'th diagonal
element is .beta..sub.mk.sup.2, the variance of the gain multiplier
.alpha..sub.mk.
Multi-Target and Multi-Scatter Models
We now generalize the model of Mutiplicative Aperture Noise to deal
with the narrow band antenna model in the presence of either
multiple scatterers or multiple emitters.
In this case we use the channel model proposed in Equation (6),
x.sub.k=H.sub.kg.sub.k+e.sub.k (10) H.sub.k=[h.sub.k1,h.sub.k2, . .
. h.sub.kQ] g.sub.k.ident.[.alpha..sub.k1,.alpha..sub.k2, . . .
.alpha..sub.kQ].sup.T g.sub.k.about.N(1,R.sub.g.sub.k.sub.g.sub.k),
where 1 is the all ones vector, and where we drop the {tilde over (
)} representation for whitened data, assuming that the data and
channels have been whitened if needed. The q index is the scatterer
index or emitter index.
From this model we obtain x.sub.k.about.N(H.sub.k1,
R.sub.e.sub.k.sub.e.sub.k+H.sub.kR.sub.g.sub.k.sub.g.sub.kH.sub.k.sup.H),
so that we can write the likelihood function as,
p(x.sub.k|H.sub.k)=det(.pi.R.sub.x.sub.k.sub.x.sub.k).sup.-1
exp(-(x.sub.k-H.sub.k1).sup.HR.sub.x.sub.k.sub.x.sub.k.sup.-1(x.sub.k-H.s-
ub.k1))
R.sub.x.sub.k.sub.x.sub.k=R.sub.e.sub.k.sub.e.sub.k+H.sub.kR.sub.g-
.sub.k.sub.g.sub.kH.sub.k.sup.H.
The mixture model in Equation (10) is very suggestive of multi-user
detection (MUD). Indeed either successive interference cancellation
(SIC) or parallel interference cancellation (PIC) are viable
techniques for separating each superimposed channel vector. The
fact that these vectors are parameterized by only a few position
and gain parameters means that it is possible to separate them
using their positions and delays as a filter. Typically the
processing will keep a set of parameters associated with a given
emitter/scatterer. These parameters can be used to generate
associated channel vectors and thus to cancel them from the current
environment using SIC or PIC. The positions parameters are then
re-estimated for each active scatterer and that enables the process
to repeat so that the cancellation can improve iteratively.
Note also that one can treat one of the "scatterers" as the static
background field. It can be canceled from the environment using the
same MUD techniques if one wants the processing to only focus on
time varying phenomenon.
Multi-Aperture Target and Phase Normalization
The Multi-target signal model in Equation (10) can be used as a
basis for a more complex and accurate target model. The receiver
can be modeled as,
.times..alpha..times. ##EQU00012## where the target channel
.alpha..sub.k is seen as a mixture of multiple related component
channels, (e.g. a human body with torso and limbs). The u.sub.k
vector is presumed to be the stationary background clutter. For
simplicity one can presume that the coefficients are Gaussian
distributed, as is the background clutter with,
g.sub.k.about.N(g.sub.k,R.sub.gg) u.sub.k.about.N(
.sub.k,R.sub.uu). The probability for the signal on and signal off
are, p(x.sub.k|ON)=N(x.sub.k-H.sub.kg.sub.k-
.sub.k,H.sub.kR.sub.ggH.sub.k.sup.H+R.sub.uu)
p(x.sub.k|OFF)=N(x.sub.k- .sub.k,R.sub.uu), from this one can
compute the posterior probability of the target present (ON) and
target not present (OFF),
.function..function..times..function..times..function..times.
##EQU00013## where p.sub.1 is the prior probability of the target
being present and p.sub.0=1-p.sub.1.
An additional modification to our model is to allow the removal of
bulk phase terms for each channel component. Thus we have,
.times..alpha..times. ##EQU00014## where c.sub.q is a unit modulus
phase equalization term. This is often necessary for smaller
wavelengths, wherein the bulk phase can change rapidly over small
changes in the target position.
One can use successive interference cancellation (SIC) or parallel
interference cancellation techniques to learn c.sub.q, simply by
rewriting,
.noteq..times..alpha..times..times..times..noteq..times..alpha..times..id-
ent. ##EQU00015##
Finally, one can use similar techniques to add and subtract channel
vectors from the collection that defines a target. The channels
that have the lowest gain parameters can be dropped and a new
channel can be added near to the canceled signal,
.apprxeq..noteq..times..alpha..times. ##EQU00016## We either pick
h.sub.p as part of the calibrated array manifold or simply as a raw
data estimate.
Position/Parameter Uncertianty
Each channel vector has a dependency on it's position and other
extraneous parameters, h.sub.kq (p.sub.q, u). It is of interest to
presume that there is some uncertainty with regards to these
parameters. For example u may contain the positions of the receiver
antennas, whilst P.sub.q may contain the uncertain
scatterer/emitter positions. We may assume or estimate priors for
these parameters, p.sub.q.about.f(p.sub.q) and u.about.g(u). If
these prior distributions are Gaussian, then we have,
p.sub.q.about.N (p.sub.0q;R.sub.p) and u.about.N (u.sub.0,R.sub.u).
Defining the real multivariate normal distribution by
(.alpha.,R).ident.det(2.pi.R).sup.-1/2exp(-1/2.alpha..sup.TR.sup.-1.alpha-
.), we can write the likelihood function as,
.function..function..times. .function..times.
.function..times..times..function..times. ##EQU00017##
Using Bayes theorem we can write out the posterior distribution
as,
.function..intg..intg..times..intg..intg..function..intg..times..intg..ti-
mes..times..intg..times..times..intg..function. ##EQU00018## While
we do not have a closed form solution for the posterior
distribution we can use a Gaussian mixture model and importance
sampling to produce a close approximation to it. These techniques
are discussed in more detail in Gaussian-Mixture-Filter,
herein.
Range Filtering
For many of the applications it is possible to know the precise
timing of the transmit waveform. This is true for example for radar
applications or for global position systems, or possibly cellular
broadcast channels. This permits the receiver to filter signals
based on their range by exploiting the time delay of arrival. This
is useful for removing unwanted multipath and other types of co
channel interference.
Consider the narrow band model where the signal is received at
multiple frequencies indexed by k. At the receiver's m'th antenna
we might see a signal of the form,
x.sub.km=h.sub.kme.sup.-2.pi.jf.sup.k.sup..tau..sup.m+i.sub.km,
where .tau..sub.m is the signal delay seen at antenna m, and where
f.sub.k is the signal frequency for index k. We perform an inverse
discrete Fourier transform, .sup.-1 and obtain,
.times..pi..times..times..times..times..times..pi..times..tau..times..pi-
..times..times..times. ##EQU00019## If its assumed that the
frequency hopping is linear f.sub.k=f.sub.0+k.DELTA., and that
h.sub.km.apprxeq.h.sub.m is dominated by a frequency independent
mode, we then have,
.times..times..pi..times..tau..times..times..pi..times..times..DELTA..fu-
nction..tau..times..pi..times..times..DELTA..function..tau.
.times..times..pi..times..function..times..DELTA..times..tau..times..func-
tion..times..DELTA..times..tau..function..DELTA..times..tau.
##EQU00020##
One can plot what the simple delay channel looks like in the time
domain. Here assuming K=128, .DELTA.=20 MHz and .tau..sub.m=5 ns,
we plot the impulsive response in FIG. 5.
As the number of frequencies K.fwdarw..infin., the time response
would approach an impulse function. For an emitter/reflector that
has an approximate known range, ie during tracking, it would be
possible to use it's round trip or total delay to filter out
unwanted multipath or interfering signals.
The required processing is illustrated in FIG. 6. The time window
represents a filtering operation in the time domain wherein only
certain delays are allowed through the filter. This is
approximately equivalent to convolving with a sinc( ) function in
the frequency domain.
The time delay filtering concept in FIG. 6 can also be used to
detect a new object or transmitter in the field of view. One sweeps
over allowed time delays and then detects any changes in the
received channel vector.
Note for other system paradigms such as multi-transmitter
geolocation, such as global positioning systems (GPS), one can do
time delay filtering simply by limiting the lags of the spreading
code correlators to a predetermined range.
Parameter Estimation and Tracking
One primary tool, used in one or more embodiments, for estimating
free parameters is Bayesian estimation. It has the feature of
providing both asymptotically optimal estimates of our parameters
and also provides post estimation error bounds. Bayesian estimation
techniques have been given a large boost by Markov Chain Monte
Carlo methods, importance sampling and particle filters, since
these techniques do not require the evaluation of intractable
integrals. Furthermore the particle techniques do not require the
computation of derivatives, or sequential optimization and are
particularly suited for multiprocessor computation, either on a
multi-core processor or an embedded graphics processing unit
(GPU).
Gaussian Mixture Particle Filters
A particularly simple and efficient particle filter technique
involves using a Gaussian or Gaussian mixture density as a prior to
generate a lot of random points. Those points are then weighted by
the conditional distribution, to then obtain a Gaussian mixture
posterior estimate. The technique can use a simple Gaussian, or to
approximate more complex priors we can use Gaussian mixtures.
While most of the prior art deal with dynamic systems and therefore
update their particles/points using the system model, we are
initially concerned with simply estimating the posterior
distribution from our set of random guesses (aka particles)
generated by our priors. Our notation below, drops the tilde for
simplicity, but the received data and channel models could, without
loss of generality, be whitened to make them interference
resistant.
Bayes theorem for computing the posterior distribution can be
written as,
.function..function..times..function..intg..function..times..function..ti-
mes. ##EQU00021## where x is the received data estimate of the
channel, and p is the parameter vector containing the location and
occurs in the estimator as the argument of h(p), the channel vector
model. If the number of received antennas is a constant M for each
channel number k, then x and h are of dimension MK. One possible
example of a likelihood function can be derived from Equation
(9),
.function..function..pi..times..times..function..times..times..times..tim-
es..ident..sigma..times..times..times..beta. ##EQU00022## As can be
seen it might be very difficult to optimize the conditional
likelihood function using ordinary techniques due to it's complex
dependency on p.
Define the real multivariate Gaussian probability density function
as, (p;R)=det(2.pi.R).sup.-1/2 exp(-1/2p.sup.TR.sup.-1p). One can
propose to model the prior and posterior probabilities as a
Gaussian mixture. The prior is assumed to adhere to,
.eta..times. .function. ##EQU00023## where
.SIGMA..sub.j.eta..sub.j=1. One can now propose the following
algorithm to update our Gaussian mixture so that it approximates
the posterior probability given in Equation (12) with
p(p)=p.sub.0(p).
Gaussian Mixture Filter Update
1. Sample N random points from the Gaussian mixture prior
p.sub.0(p), p.sub.n.
2. Compute the conditional likelihood function over all the random
points, .lamda..sub.n=p(x|p.sub.n).
3. Normalize the likelihoods to sum to one,
.lamda..fwdarw..lamda..lamda. ##EQU00024##
4. Resample a subset of Q<N points from the original N, where
the probability of choosing p.sub.n equals .lamda..sub.n. Call this
set of renumbered points p.sub.v(q) for subindex q and the
associated original index, n=v(q).
5. Use the EM algorithm, and the points p.sub.v(q) to learn a new
set of n'.sub.j, m'.sub.j and R'.sub.j to fit the set of sampled
p.sub.v(q) points, q=1 . . . Q.
We now delve into some of the details of the algorithm.
Sampling the Mixture and Resampling
For a given mixture probability of the form,
.eta..times. ##EQU00025## we can generate N pseudo-random samples
from any such distribution for which we can sample directly from
f.sub.j(p). For each n we randomly select what category j, p.sub.n
will belong to, based on the category probabilities .eta..sub.j.
Selecting the category can be done in approximately o(log(J))
comparisons by using a binary search into the cumulative sum array
of .eta..sub.j, namely
c.sub.j.ident..SIGMA..sub.q=1.sup.j.eta..sub.j. When a uniform
random variate u is chosen, we need merely find j such that,
c.sub.j-1.ltoreq.u<c.sub.j, using binary search, and
initializing c.sub.0=0.
Once the category is chosen for a given n, we generate a random
sample using the pdf f.sub.j(p). For example when
f.sub.j(p)=(p-m.sub.j;R.sub.j), we can sample from any given
multivariate Gaussian distribution by setting
p.sub.m=m.sub.jC.sub.j.sup.Te.sub.n, where e.sub.n is the n'th,
generated unit variance Gaussian random vector with independent
components, each component selected according to any number of
standard proceedures for generating Gaussian samples e.g. the
ziggurat method. Here the C.sub.j matrix is the upper Cholesky
factor of the positive definite covariance matrix R.sub.j. These
computations are highly amenable to parallel processing as well
After this we choose to sub-sample our original large sample of N
points in proportion to the likihood
.lamda..function..function. ##EQU00026## Again we can compute a
cumulative sum of the likelihoods
v.sub.n.ident..SIGMA..sub.q=1.sup.n.lamda..sub.n, followed by a
binary search to find a n where a given, newly generated random
uniform variable u.sub.q satisfies, v.sub.n-1.ltoreq.u<v.sub.n.
This is our n=v(q).
Let us compute the probability distribution for p.sub.v(q). We
have,
.function..function..function..times..function..function..times..lamda..f-
unction..times..noteq..function..function..apprxeq..function..times..times-
..intg..function..times..times..function..times..intg..function..times..ti-
mes..function. ##EQU00027## We have
.SIGMA..sub.q.noteq.np(x|p.sub.q)+p(x|p).apprxeq.N
.intg..sub.pp(x|p)p.sub.0(p)dp, by the law of large numbers, i.e.
the sample mean approaches the expected mean as the number of
independent samples N, approaches infinity, and p is also chosen
independently. This shows that the conditional distribution of
p.sub.v(q) asymptotically approaches our desired posterior
distribution. In what follows we label the resampled points as
p.sub.v(q).fwdarw.p.sub.q.
Expectation Maximization Algorithm
The Expectation Maximization (EM) algorithm can be used to learn
the free parameters of a Gaussian mixture model, from collected
data. In Samling the Mixture and Resampling it was showed how to
approximately sample from the posterior distribution given an
arbitrary likelyhood function. Thus this is the final step needed
to update the parameters associated with the posterior distribution
p(p|x).
The Gaussian mixture model for the posterior distribution will take
the form,
.pi..times. .function. ##EQU00028##
The initial values for .pi..sub.j, m.sub.j, and R.sub.j are chosen
to be the same as the prior. The EM algorithm consists first
performing the expectation of the logarithm of the conditional
distribution p.sub.1(p|m,R), which is a lower bound to the original
function followed by a maximization of this lower bound over the
unknown parameters m.sub.j and R.sub.j, to find the next iterate.
This can be summarized as provided below.
Expectation Step
Compute the posterior likelihood of a sampled point p.sub.q being
in category j,
.gamma..pi..times..times. .function..pi..times.
.function..times..ident..gamma. ##EQU00029## The parameter N.sub.j
can be thought of as the expected number of points associated to
category j. Note .SIGMA..sub.jN.sub.j=Q.
Maximization Step
The updated mixture parameters are given by,
.times..gamma..times..times..times..gamma..times..times..pi.
##EQU00030## The Expectation and Maximization steps are performed
repeatedly until the log-likelihood function increase is below a
predefined threshhold. The log likelihood function is defined
as,
.function..pi..times. .function. ##EQU00031##
The form of EM algorithm suggests a modification that potentially
allows one to skip the resampling component of Sampling the Mixture
and Resampling. Using the original sampling of p.sub.n drawn from
the prior distribution, modify (13) to read,
.gamma..pi..times. .function..times..function..pi..times.
.function..times..function. ##EQU00032## This causes a given
p.sub.n to be directly weighted by it's likelihood function.
Geolocation and Tracking
We can use our Bayesian particle filter techniques to facilitate
both geolocation and tracking. For a single target objective we can
use Equation (9) as an objective function to perform an initial
maximization over position to find p. The position parameter enters
the equation through the channel model h(p).
Dynamic Models
One simple way to approach dynamic models, is to simply re-estimate
changing parameters and allow them to deviate from their current
values at each look at the environment. This can be accommodated by
our Bayesian approach, simply by bounding the variance of the
parameter of interest from below. So for example if .rho.(p|x) is
the posterior distribution of position parameter p given the
observed data x, after many observations we might have,
.rho..function..rho..function..times. ##EQU00033## We then can
exponentially devalue stale information by updating in the
log-domain, log(p.sub.N(p|x))=log(.rho.(p|x.sub.N))+.lamda.
log((p.sub.N-1(p|x))+c, where 0<.lamda.<1 and c is a
normalization constant so that p.sub.N(p|x) integrates to 1. The
exponential average puts a lower bound on the variance, which would
tend towards zero if all the prior information was weighted
equally. This has the effect of allowing a search over the
parameter space in the vicinity of the current optimal p, with
enough variance so that p can change according to system dynamics.
This model works well for slow moving reflectors/emitters.
Another approach for system dynamics is to assume that p adheres to
a first order discrete constant velocity model, ie
p.sub.N=p.sub.N-1+v. The model then seeks to learn v and to update
p.sub.N based on the observed statistics. Again we can limit the
variance for estimating v so that we can allow the acceleration to
change slowly over time.
Every emitter q in the environment, will have an associated set of
parameters that will need to be tracked and maintained, aka
p.sub.q. Moreover for some systems, e.g. detecting hand gestures,
there may be multiple emitters/reflectors associated with the same
target, ie a hand arm and human body. Emitters and reflectors can
typically be associated by position over time relative to one
another. We can use machine learning techniques to associate such
objects and to classify their trajectories using channel vectors as
the training data. A Bayesian belief network is a good model for
this.
Detecting Environment Changes
During normal operation of any proposed sensor device, the
environment may be quasi-stationary. Many applications require
detection of moving or living objects and thus are improved by
removing stationary objects. To this end we can consider a
hypothesis channel model of the form, x.sub.off=u+e.sub.0
x.sub.on=u+h+e.sub.1, where the channel vectors are stacked over
all spread frequencies/look indices and have length MK elements,
where u is the stationary receive vector due to either stationary
emitters or stationary scatterers, e.sub.1 is the noise vector for
the new signal present case, e.sub.0 is the signal off noise, and h
is the channel vector for a new or moving emitter.
Because we continuously monitor the array we can learn an
approximate prior distribution for u and h wherein u.about.N (
,R.sub.uu) and h.about.N (h,R.sub.hh). The statistics of u can be
updated over each quiescent period. Suppose that we receive new
data x. It is apparent that the conditional likelihood of x given
OFF or ON, can be written as,
p(x,u|ON)=N(x-u-h,R.sub.hh+R.sub.ee)N(u- ,R.sub.uu)
p(x,u|OFF)=N(x-u,R.sub.ee)N(u- ,R.sub.uu). After integrating out u
one obtains, p(x|ON)=N(x-h- ,R.sub.uu+R.sub.hh+R.sub.ee)
p(x|OFF)=N(x- ,R.sub.uu+R.sub.ee) From this we can see that the
posterior probability of x adhering to the ON model is,
.function..function..times..function..times..function..times.
##EQU00034## where p.sub.1 is the prior probability that we see a
new signal arriving and p.sub.0=1-p.sub.1. In general the
distribution for h will be a Gaussian mixture, conditioned on the
current set of captured targets and the possibility of a new
acquisition. Thus we may choose to further condition the
distribution on a category k, that defines the existing tracked
targets plus an additional zero mean distribution for an unknown
new target.
We also need to compute the posterior distribution for u so that we
can update it's statistics. We have, p(u|x,ON)=N( .sub.1,R.sub.11),
.sub.1= +R.sub.uu(R.sub.uu+R.sub.hh+R.sub.ee).sup.-1(x- -h)
R.sub.11=R.sub.uu-R.sub.uu(R.sub.uu+R.sub.hh+R.sub.ee).sup.-1R.sub.uu;
p(u|x,OFF)=N( .sub.0,R.sub.00), .sub.0=
+R.sub.uu(R.sub.uu+R.sub.ee).sup.-1(x- )
R.sub.00=R.sub.uu-R.sub.uu(R.sub.uu+R.sub.ee).sup.-1R.sub.uu. The
preferred modality is to update only when the OFF state is
detected. Note also the similarity of these updates to the Kalman
filter. The latter has a close relationship to Bayesian inference,
when the distributions in question are normally distributed.
Shape and Gesture Identification
Suppose that a target has a collection of emitters/scatterers
associated with it h.sub.q, q.di-elect cons. where e is a target
index and is a set of scatterers associated with a given target.
The shape identification problem is to identify the set of channel
vectors h.sub.q that correspond to a shape category c.sub.e. The
identification problem also often has a time dependency, so that
one typically makes an identification over a time period t.di-elect
cons.[r-T,.tau.]. An example of this would be say a human hand
making a right circular motion. To identify this, we require data
over a period of time long enough to verify that it is a circular
motion in the right direction.
Mapping a set of channel vectors to a shape is a good candidate for
machine learning and using a Bayesian belief network. The training
process in learning and identifying targets is part of the
calibration process for the array. Each target will have very
different spatial signatures, especially after the target is
observed over time. The category associated with a target's shape,
the target's orientation as well as it's position are all candidate
parameters for learning within the p vector of the overall
parameter estimation problem.
While the data set is large, the identification problem can often
be broken down into a set of smaller more reasonable steps. For
example we can consider the problem of detecting a human hand wave.
If we model body parts as cylinders, we need to find a cylinder
"object" making changes in the elevation coordinate, or more
precisely a cylinder that is rotating about an axis in the x-y
plane. The gesture detection problem is roughly outlined in FIG.
7.
The initial state occurs when the target of interest appears
stationary or has not yet entered the field of view. Once motion is
detected via the techniques in Detecting Environment Changes we
then use Bayesian inference to detect the object location, it's
orientation and object shape using the Bayesian inference
techniques of Gaussian Mixture Pareticle Filters and Geolocation
and Tracking.
For the gesture detection we pay attention to it's orientation and
it's elevation coordinate. After two changes in elevation
coordinate we can presume that the object is waving. Not shown in
FIG. 7 are all the possible transitions to the stationary state,
which can happen after a time-out period in any give state. When
the hand wave gesture is detected, that can be reported to the
overall system controller.
Many other gestures and motion detection schemes can be broken up
into discrete steps, whose transitions can be appropriately
detected by Bayesian inference and machine learning in general. It
is important for us to work directly with the channel vectors,
since this gives us the optimal theoretical estimation and
detection performance.
Using HOSVD in Machine Learning
Creating a neural network, Bayesian network or other machine
learning structure, allows one to learn an arbitrary nonlinear
function through a controlled data structure. An example of such a
structure is provided in FIG. 8.
In this directed graph each node represents a computation. For a
Bayesian network the graph represents conditional dependency and
each edge is weighted by the conditional probability, and the nodes
represent a sum of products whose output is the posterior
probability, ie for node k,
.di-elect cons..function..times. ##EQU00035## The initial
probabilities and conditionals are possibly dependent on the
received input data and par(k) is the set of nodes which are
parents of node k in the directed graph.
For a neural network the weights on each edge multiply the values
at the nodes are summed and then passed through a nonlinear
function to obtain the node output,
.di-elect cons..function..times. ##EQU00036## The machine learning
literature contains many techniques for learning the weights and
the structure of these networks in an attempt to match an arbitrary
nonlinear function.
We introduce here a new approach that exploits the structure
compression properties of higher order singular value
decomposition. In FIG. 9 we show that the HOSVD network is
essentially a one layer network that learns a set of
eigen-functions (the f.sub.kj) and some multipliers to approximate
the nonlinear function of interest. The functions are multiplied
and then summed to obtain the final function value via,
.function..apprxeq..sigma..times. ##EQU00037## The number of
branches required to approximate most functions is significantly
smaller than the requirement to tabulate the multivariable
function. For example, when we look at HOSVD decompositions for
Green function kernels in Basis Functions and HOSVD and Equation
(31) we find only a tiny fraction of the singular values are
significant as FIG. (11) shows.
Thus there may be significant computational savings in using this
flattened structure, once the individual eigenfunctions are
approximated by splines. This structure provides an innovative way
of implementing many of the machine learning networks and for
approximating the many functions used in partial differential
equation boundary value problems, such as the Helmholtz wave
equation for electromagnetic propagation.
Channel Modeling
A key component of this invention is the ability to successfully
model the electro-magnetic channel seen by our receivers.
Particularly challenging is the problem of electro-magnetic
scattering. We propose a series of approximations to the scattering
problem that exploit known geometries and some novel expansions of
the electromagnetic reproducing kernel.
Our goal is to create a series of basis functions that can be used
to calibrate our receiver array in the presence of various
scattering objects. These functions can also be used to interpolate
our scatterers over differing range and angle of arrival
values.
Electro-Magnetic Scattering
For our calibration model, we desire to pick a set of basis
functions that match the electromagnetic environment of the
receiver. This is especially important for nearby scatterers for
radar applications. There are some applications that require us to
model a sizeable reflector within a few feet of the transmitter.
For example for security purposes we might want to differentiate
between an animal or human reflector, or we may want to detect hand
or arm gestures or detect the presence of a forklift entering a
warehouse etc.
For this reason we need to extend the emitter model beyond that of
a point source. If we model the scatterer as a parameterized
surface that admits a local orthogonal coordinate system, we can
exploit a type of Green's theorem to model an electromagnetic field
as a function of the field parameters on the distant closed
surface. Indeed it is well known that the wave equation has a
reproducing kernel, similar to Cauchy's theorem in complex
analysis.
In the geometric Clifford algebra Cl.sub.0,3, the microscopic
Maxwell's equations can be written as,
.times..differential..times..times..epsilon..times..times..rho.
##EQU00038## where c is the velocity of light, .differential..sub.t
is the partial derivative operator with respect to time,
D.ident.e.sub.1.differential..sub.x+e.sub.2.differential..sub.y+e.sub.3.d-
ifferential..sub.z is the Dirac operator, F.ident.E+icB is the
electromagnetic field vector (or Faraday vector), .rho. is the
charge density and j is the current density vector. A vector is
expressed by v.ident..SIGMA..sub.k=1.sup.3v.sup.ke.sub.k, and the
coordinate vectors e.sub.k have an associative product that has the
properties e.sub.qe.sub.m=-e.sub.me.sub.q for q.noteq.m and
e.sub.q.sup.2=1. The pseudoscaler i.ident.e.sub.1e.sub.2e.sub.3
commutes with the e.sub.k and has the structural property
i.sup.2=-1.
A product of vectors in the geometric algebra includes both the
inner product and the cross product via the formula,
ab=ab+ia.times.b. Because of this it is possible to express
Maxwell's equations in the concise expression 14. In a region of no
sources, and assuming the time dependency is harmonic, of the form
exp(2.pi.jft), we see that the electromagnetic field vector
satisfies the Helmholtz equation, DF=-jkF, where the wavenumber
.times..times..pi..times. ##EQU00039## Note that the j used here is
a square root of -1 for tracking phase changes for a complex phasor
at frequency f, which is different from the pseudoscalar i. An
important property of the Dirac operator within a Clifford algebra
is the fact that,
D.sup.2=.differential..sub.x.sup.2+.differential..sub.y.sup.2+.diff-
erential..sub.z.sup.2 is the scalar Laplace operator of the
standard wave equation.
The Cauchy integral formula for our Geometric algebra can be
written in general by,
.intg..differential..times..intg..times..times. ##EQU00040## where
dx*.ident.dydze.sub.1-dxdze.sub.2+dxdye.sub.3=.SIGMA..sub.kdx*.sub.ke.sub-
.k, and * is the Hodge star operator, R is a volume in three space
and .differential.R is the boundary of the volume. The vector form
dx* will be normal to any parameterized surface. In (15) if we
substitute h.fwdarw.b and g.fwdarw.Da, we get,
.intg..differential..times..intg..times..times..times. ##EQU00041##
Similarly if we substitute h.fwdarw.bD and g.fwdarw.a we get,
.intg..differential..times..intg..times..times..times. ##EQU00042##
Subtracting these two equations yields a second order Green's
formula,
.intg..differential..times..times..intg..times..times..times..times.
##EQU00043##
Note that it is possible to extend both the Cauchy equation and
Green's formula to the time harmonic case, addressing the Helmholtz
equation. This is simply a matter of adding and subtracting the
eigen value, or simply letting them cancel out for the Green's
function. We can thus write,
.intg..function..times..function..times..times..intg..differential..times-
. ##EQU00044##
.intg..differential..function..times..times..times..times..intg..times..t-
imes..function..times..times. ##EQU00044.2##
Bicomplex Algebra
The addition of another commutative square root of negative one, j,
for the purpose of tracking harmonic phase, causes the underlying
geometric algebra to split into 2 isomorphic components. One
component corresponds to outgoing waves with a wavenumber k, and
the other to incoming waves with wavenumber -k. It also forces
plane waves to decompose into right and left circular
polarization.
The two algebras coincide with coefficients of either
.sub.+=1/2(1+ij), and .sub.-=1/2(1-ij). The algebra isomorphism is
better understood by noting, .sub.+.sub.+=.sub.+
.sub.-.sub.-=.sub.- .sub.+.sub.-=0 .sub.++.sub.-=1
i.sub.--i.sub.+=j .sub.-e.sup.iwt+.sub.+e.sup.-iwt=e.sup.jwt. Thus,
we can write the Clifford electromagnetic field as,
.function..times..times..omega..times.
.times..times..times..times..omega..times.
.times..times..times..omega..times.
.times..function..times..times..omega..times.
.times..function..times..times..omega..times. ##EQU00045## where
the bar here is conjugation with respect to j. The j is replaced by
i and it's conjugate in the split algebra. This isomorphism allows
one to drop the additional j and implement the algebra in our
calculations by tracking Clifford numbers for both
.times..pi..times..omega..times..times..times..pi..times..omega.
##EQU00046## One can recover the original representation by
replacing, .sub.+.fwdarw.1/2(1+ij) .sub.-.fwdarw.1/2(1-ij),
resulting in,
.function..times..times..omega..times..times..function..times..times..ome-
ga..times..times..times..omega..times..times..times..function..times..time-
s..omega..times..times..times..omega..times..times..times..omega..times.
##EQU00047## where, F.sub.r=1/2(F(k,r)+F(k,r))
F.sub.i=1/2(F(k,r)-F(k,r)).
Boundary Conditions
For the scattering problem it is important to understand
electromagnetic boundary conditions. We presume linear media so
that Equation (14) holds. In Equation (16) set h=1 and let g=F, the
electro-magnetic field vector, so that (D+jk)
.times..epsilon..times..times..rho. ##EQU00048##
Let the volume of integration be an infinitely thin box that spans
the boundary of two regions as shown in FIG. 10. Let A be the face
of the box in free space and A' the region on the other side of the
boundary surface. We thus have,
.times..epsilon..times..times..rho..intg..intg..differential..times..intg-
..differential..function. ##EQU00049## where .rho..sub.s is the
surface charge density and j.sub.s is the surface current density
and n is the outward facing normal from surface A in region 1 to
the left, and is negative to the normal in region 2 to the right of
the boundary. dS is the scalar area differential. Since F is
piecewise continuous and the volume of R.fwdarw.0, we also have,
.intg..sub.RjkF.fwdarw.0. The sides of the box shrink to 0 rapidly
and don't contribute to the area integral. The A region is presumed
small enough so that F is constant there.
In one fell swoop, this proves the electromagnetic boundary
conditions,
.function..epsilon..times..rho..times..mu..times. ##EQU00050## This
shows in particular that, n.times.(E.sub.1-E.sub.2)=0
n(B.sub.1-B.sub.2)=0 The tangential component of the electric field
is continuous, as is the normal component of the magnetic field. We
can project the Faraday vector on to it's tangential E component
and normal B component via the surface projection operator,
Q.sub.-F.ident.1/2(F+n{tilde over (F)}n), where a is the principle
involution defined by {tilde over (e)}.sub.k=-e.sub.k for k=1 . . .
3, and ab=a{tilde over (b)}. Thus the continuity relations can be
written as, Q.sub.-(F.sub.1-F.sub.2)=0. If you define the
complementary projection operator onto the normal component of E
by, Q.sub.+F=1/2(F-n{tilde over (F)}n), then we have on the
surface,
.times..epsilon..times..rho..times..times..mu..times..times..function..ep-
silon..times..rho..times..mu..times. ##EQU00051##
For a perfect conductor the surface current paravector
.epsilon..times..rho..times..mu..times. ##EQU00052## can be found
from the incident Faraday vector F.sub.0 by, K=2nF.sub.0, where is
the real part of the Clifford number a with respect to the
pseudoscalar i=e.sub.1e.sub.2e.sub.3. Here the fields are evaluated
at the surface of the conductor. We can also write down the
reflected field from a perfect conductor as, F.sub.r(x)=-n{tilde
over (F)}.sub.0(n{tilde over (x)}n)n, Also of interest is the
reversion involution defined by
(ab).sup..dagger.=b.sup..dagger.a.sup..dagger. and
e.sub.k.sup..dagger.=e.sub.k. Note that =-i, and i.sup..dagger.=-i.
Also a=1/2(a+a.sup..dagger.). For non-deal conductors we will
generalize this rule somewhat via an unknown reflection coefficient
.alpha., so that, F.sub.r(x)=-.alpha.n{tilde over
(F)}.sub.0(n{tilde over (x)}n)n (18) The reflection coefficient
.alpha. is allowed to be complex in j but not in i. These formulas
require that the coordinate system is chosen so that the reflecting
conducting plane passes through the origin. Otherwise we need to
replace n{tilde over (x)}n with n(x-x.sub.0){tilde over (
)}n+x.sub.0, for some x.sub.0 on the reflecting surface.
This model yields the following surface current paravector (vector
plus a scalar), K=2nF.sub.0.alpha..sup..dagger. (19) where F.sub.0
is evaluated at the reflector surface. For monochromatic waves it
suffices to evaluate the Faraday vector over a closed surface to
determine it's radiated pattern throughout space. This makes
Equation (19) particularly useful to determine the electromagnetic
field, reflected off of arbitrary surfaces.
It can be shown that the reflected electromagnetic field vector in
Equation (18) satisfies Diracs equation,
.times..differential..times..times..times..omega..times.
##EQU00053## since the reflection formula is in fact just an
improper Lorentz transformation. We can use Equation (19) as an
approximation for non-ideal conductors, wherein the complex
reflection coefficient .alpha. is a property of the underlying
material. Since .alpha. can be learned during calibration and
estimated from the received data, changes in .alpha. can be
attributed to changes in the underlying dielectric material of
different objects. This has applications to medical imaging, where
the dielectric may change as a result of sugar levels, or tumors
and so forth.
The Reproducing Kernel
Consider the Helmholtz form of the Cauchy equation in (16). Suppose
we can find a function h(x) in the geometric algebra that has the
property that applying the right Helmholtz operator yields the
Dirac delta function, h(x-x.sub.0)(D-jk)=.delta.(x-x.sub.0).
Furthermore consider a Faraday vector g(x) satisfying the Helmholtz
version of Maxwell's equation,
.times..function..times..epsilon..times..times..rho. ##EQU00054##
Substitution into (16) yields,
.intg..delta..function..times..function..times..times..epsilon..times..ti-
mes..rho..function..function..times..intg..differential..function..times..-
times..function. ##EQU00055## Assuming that x.sub.0 lies in the
spatial region R (possibly unbounded assuming decaying radiation
conditions), we have,
.function..intg..differential..function..times..times..function..intg..fu-
nction..times..times..epsilon..times..times..rho..function..function..time-
s. ##EQU00056## For the case of no source charges or currents in
the field of interest, we have an analogy to Cauchy's equation in
complex analysis,
.function..intg..differential..function..times..times..function.
##EQU00057## The reproducing kernel for geometric algebra can be
written as,
.function..times..times..pi..times..ident. ##EQU00058## This can be
written as,
.function..times..pi..times..times..times. ##EQU00059## .ident.
##EQU00059.2##
For the Green's theorem in (17), it is helpful to first define the
Electromagnetic paravector potential as, a=cA-.PHI., where A is the
electro-magnetic vector potential and where .PHI. is the scalar
potential. We also assume the Lorenz gauge which forces
.differential..differential..times. ##EQU00060## where .sub.S is
the scalar part of a Clifford number. This means,
.times..differential..PHI..differential. ##EQU00061## Maxwell's
equations now give us,
.differential..differential..times. ##EQU00062## For the time
harmonic case we get the Helmholtz version of this as,
.times..times..times..times..epsilon..times..times..rho.
##EQU00063## In Equation (17) we now use a as the electromagnetic
paravector potential and b as a Green's function satisfying
(D.sup.2+k.sup.2)b (x-x.sub.0)=.delta.(x-x.sub.0). This yields,
.intg..differential..function..times..times..times..function..function..t-
imes..times..function.
.times..intg..function..times..epsilon..times..times..rho..function..func-
tion..times..function. ##EQU00064##
It is customary to choose a Green's function that is constant or
even zero on .differential.R, however our typical use case
evaluates the paravector potential at an arbitrary point x.sub.0 in
space, requiring our surface to vary with the choice of receiver
location x.sub.0. For this reason it suffices to simply choose the
standard Green's function for the sphere,
.function..times..pi..times. ##EQU00065## from which we can
obtain,
.times..times..ident..function..times..times..pi..times..times..times.
##EQU00066## In a region with no sources we can therefore write the
reproducing kernel for the paravector potential as,
.function..intg..differential..times..times..times..times..function..func-
tion..times..times..function. ##EQU00067## Note also that if the
region R is unbounded and that we assume a radiation decay
condition that makes the surface integrals vanish, we obtain the
usual time harmonic solution for the retarded potential,
.intg..function..times..epsilon..times..times..rho..function..function..t-
imes..function. ##EQU00068## Also note that the scalar potential
can, in particular can be written as,
.PHI..function..intg..differential..times..times..times..times..function.-
.function..times..times..function. ##EQU00069##
Reproducing Kernel Channel Model and Surface Inference
The reproducing kernel formulas suggest a possible numerical model
for the reflected waveform. We pick a hypothesized surface R
centered on the reflector position p, and with boundary
.differential.R. Given a model for the electromagnetic field from
the transmitter, we can compute the reflected field from Equation
(18), based on the local normal to the surface .differential.R,
from this we can compute the field at the remote receiver antenna
using one of the reproducing kernel formulas in Equation (20) and
Equation (24).
We summarize this algorithm as provided below.
Reproducing Kernel Algorithm
1. Hypothesize a position p for the center of the emitter/reflector
region R.
2. Compute a set of Faraday vectors for the reflected field
.alpha..sub.qF.sub.r(x.sub.q)=-.alpha..sub.qn.sub.q{tilde over
(F)}.sub.0(p+r.sub.q)n.sub.q, where F.sub.0(x) is the transmitted
Faraday vector seen at position x, n.sub.q is the unit normal
vector seen at position x.sub.q=p+r.sub.q.di-elect
cons..differential.R.
3. Compute the Faraday vector at antenna position y using an
approximation to the reproducing kernel formula,
.function..apprxeq..function..times..times..alpha..times..times.
##EQU00070## where A.sub.q is the area associated with the sampled
surface point x.sub.q.
4. Compute per antenna voltage h.sub.m.apprxeq.b.sub.m E(y.sub.m),
where b.sub.m is the appropriate vector effective length for
receiver antenna m, and E(y.sub.m) is the electric field component
of g(y.sub.m). Computations can be reduced by bringing the E field
projection and the inner product into the sum/integral, ie
.apprxeq..alpha..times..times..function..times..times. ##EQU00071##
where u.sub.RV is the real vector part of u. Of course more
sophisticated quadrature formulas can be used to obtain the
integral approximation.
5. Use the computed channel vector h.sub.m as part of the Gaussian
Mixture Filter Update algorithm from Gaussian Mixture Filter Update
for estimating the reflector/emitter position p.
Another algorithm suggests itself from the form of the Reproducing
Kernel Algorithm. It is theoretically possible to compute the
contribution of all the nodes to the m'th antenna via,
.sub.mq.ident.A.sub.qb.sub.mh(x.sub.q-y.sub.m)n.sub.qF.sub.r(x.sub.q).sub-
.RV, in advance after calibration. In this mode we would compute
.sub.mq over a dense grid x.sub.q, indexed by q. Thus with the
.sub.mq pre computed we can write an approximation for the channel
vector as,
.di-elect cons..differential. .times..alpha. ##EQU00072## In this
model we might consider the reflection coefficients .alpha..sub.q
and even the boundary set .differential.R as free parameters, not
to mention the implicit dependency on the surface normal unit
vectors n.sub.q. By making these parameters invariant or part of
some constrained manifold we might actually be able to infer the
boundary from the collection of the received data over time.
Allowing the boundary to change permits us to possibly detect human
gestures, breathing or more complex objects in the field. The free
parameters can be added to the global invariant parameter vector p
and optimized using Gaussian Mixture Filtering.
This suggests the algorithm as provided below.
Surface Inference Algorithm
1. Hypothesize a boundary region .differential.R, containing a
center as part of the parameter vector p.
2. Identify all nodes q.di-elect cons..differential.R and compute
or update the unit normal vectors n.sub.q if necessary, based on
the shape of .differential.R.
3. Compute h.sub.m=.SIGMA..sub.q:x.sub.q.sub..di-elect
cons..differential.R .sub.mq.alpha..sub.q and use it to compute the
likelihood for use in (9) and the algorithm in the Gaussian Mixture
Filter.
4. Use the Gaussian Mixture Filter algorithm to update the mixing
gains .alpha..sub.q in a locally constrained manner (ie small
variance), and use nearby q to determine if the bondary
.differential.R should be updated if it achieves a better fit.
5. Test convergence if no improvement can be made in the
likelihood, if not return to step 1.
Factoring Common Phase Center
For higher frequencies (larger k), spherical harmonics produce
phase terms of the form e.sup.-jkr, where r is the total distance
traveled for the transmitted waveform (about twice the distance to
the reflector). The
##EQU00073## wave solution is always asymptotically true in the far
field. Because r is usually much larger than the inter-element
spacing of the array, in order to avoid floating point overflows or
underflows it helps to factor out the phase response due to the
center of the array.
Suppose the reflector is at position y and the center of the array
is at x.sub.0, which we also assume is the transmitter position for
simplicity. Assume also that x.sub.m=x.sub.0+d.sub.m is the
position of receive antenna m. The total distance traveled to
antenna m is,
r.sub.m=y.sub.0+.parallel.x.sub.0-y+d.sub.m.parallel.. The total
distance to x.sub.0 is,
.times..ident..times. ##EQU00074## If we factor out the phase
center we have the relative phase term,
e.sub.rel.ident.e.sup.-jk(r.sup.m.sup.-r.sup.0). But now we can
write,
.times..times. ##EQU00075##
The last expression avoids the loss of precision inherent in
subtracting two large numbers r.sub.m-r.sub.0 in the phase
response. These terms are inherent in all models of
electro-magnetic propagation, and it is often beneficial to work
with the relative phase. Indeed one can smooth the array response
by tossing away the bulk phase term e.sup.-jkr.sup.0 in the channel
models via relative gain normalization.
Analytic Scattering Models
We now turn our attention to techniques which approximate the
electro-magnetic field with analytic solutions using known
functions or series of functions. We seek to obtain sets of basis
functions that can be used to fit the scattering model at hand,
usually as part of the calibration process.
We first look at some useful coordinate systems, whose level curves
represent objects of interest. We also derive their Dirac
operators.
Rectangular Coordinate System
.times..differential..differential..differential..differential..times..di-
fferential..differential. ##EQU00076## This coordinate system is
good for representing rectangular surfaces and boxes.
Cylindrical Coordinate System
.rho..function..theta..times..rho..function..theta..times..rho..times..th-
eta..times..times..rho..differential..differential..rho..times..theta..tim-
es..differential..differential..theta..times..differential..differential..-
times..rho..function..theta..times..function..theta..times..times..times..-
theta..function..theta..rho..times..function..theta..rho..times.
##EQU00077## This coordinate system is good for representing
cylindrical surfaces.
Spherical Coordinate System
.times..function..theta..times..function..PHI..times..times..function..th-
eta..times..function..PHI..times..times..function..theta..times..times..th-
eta..times..PHI..times..differential..differential..times..theta..differen-
tial..differential..theta..times..PHI..differential..differential..PHI..ti-
mes..function..theta..times..function..PHI..times..function..theta..times.-
.function..PHI..times..function..theta..times..times..times..theta..functi-
on..theta..times..function..PHI..times..function..theta..times..function..-
PHI..times..function..theta..times..times..times..PHI..function..PHI..time-
s..function..theta..times..function..PHI..times..function..theta..times.
##EQU00078## This coordinate system is good for representing
spheres.
Prolate Spheroidal Coordinate System
.times..sigma..times..tau..times..function..PHI..times..times..sigma..tim-
es..tau..times..function..PHI..times..times..sigma..times..tau..sigma..tim-
es..times..times..tau..times..times..times..PHI..times..times..sigma..diff-
erential..differential..sigma..times..tau..differential..differential..tau-
..times..PHI..differential..differential..PHI..times..times..sigma..functi-
on..times..times..times..function..times..times..times..function..function-
..times..times..times..times..times..tau..function..times..times..times..f-
unction..times..times..times..function..function..times..times..times..tim-
es..times..PHI..function..PHI..times..sigma..times..tau..times..function..-
PHI..times..sigma..times..tau..times..ident..times..ident.
##EQU00079## Prolate spheroidal coordinate systems are good for
modeling an ellipse that has been rotated about the z-axis.
Confocal Ellipsoidal Coordinates
This confocal ellipsoidal coordinate system has level surfaces
corresponding to ellipsoids when .xi..sub.1 is constant,
.xi..xi..xi. ##EQU00080## hyperboloids of one sheet when .xi..sub.2
is constant,
.xi..xi..xi. ##EQU00081## and hyperboloids of two sheets when
.xi..sub.3 is constant
.xi..xi..xi. ##EQU00082## where we additionally require,
.xi..sub.1<c.sup.2<.xi..sub.2<b.sup.2<.xi..sub.3<.alpha..s-
up.2. Coordinates can be derived as,
.xi..times..xi..times..xi..times..times..xi..times..xi..times..xi..times.-
.times..xi..times..xi..times..xi..times. ##EQU00083## The Dirichlet
operators for the confocal ellipsoidal coordinates must be derived
implicitly since (.xi..sub.1,.xi..sub.2,.xi..sub.3) can only be
found from the roots of the third order polynomial equation in
z,
##EQU00084##
Cauchy-Kovalevskaya Extension Theorem
Similar to reproducing kernel concept, in the Reproducing Kernel
portion of this document. it is possible to generate an analytic
solution to the wave equation, given it's response on some closed
surface. Suppose that the Dirac operator in a given coordinate
system can be written as,
.times..psi..differential..differential..psi..times..psi..differential..d-
ifferential..psi..times..psi..differential..differential..psi.
##EQU00085## for a coordinate parameterization
(.psi..sub.1,.psi..sub.2, .psi..sub.3). Suppose we know that given
a fixed .psi..sub.3=0, the coordinate system forms a closed surface
(e.g. for spherical coordinates if .psi..sub.3=r-1, then the closed
surface is a sphere).
Define,
.times..psi..differential..differential..psi..times..psi..differential..d-
ifferential..psi. ##EQU00086## Furthermore suppose we have a
differential operator .mu.(.psi..sub.3) that depends on .psi..sub.3
and let .PHI.(.psi..sub.1,.psi..sub.2) be the value of F restricted
to the surface .psi..sub.3-0. We will make the assumption that,
(D+jk).mu.(.psi..sub.3).PHI.(.psi..sub.1,.psi..sub.2)=0 (26) so
that
F(.psi..sub.1,.psi..sub.2,.psi..sub.3)=.mu.(.psi..sub.3).PHI.(.psi..sub.1-
,.psi..sub.2), satisfies the Helmholtz equation. If we require the
operator .mu.(.psi..sub.3) to be valid for any surface function
.PHI.(.psi..sub.1,.psi..sub.2), then Equation (26) requires the
following functional equation,
.times..psi..differential..differential..psi..times..mu..function..psi..t-
imes..psi..times..mu.'.psi..times..mu..function..psi..mu.'.psi..times..psi-
..times..times..mu..function..psi. ##EQU00087## This yields the
operator solution to .mu.(.psi..sub.3),
.mu.(.psi..sub.3)=exp(-D.psi..sub.3).sup.-1D.sub.12.psi..sub.3).
Thus the solution to the Helmholtz equation becomes,
F(.psi..sub.1,.psi..sub.2,.psi..sub.3)=exp(-(D.psi..sub.3).sup.-1D.sub.12-
.psi..sub.3).PHI.(.psi..sub.1,.psi..sub.2). As long as
.PHI.(.psi..sub.1,.psi..sub.2) is C.sup..infin., we can expand the
operator function exp(-(D.psi..sub.3).sup.-1D.sub.12.psi..sub.3) in
power series of functions and derivatives with respect to
.psi..sub.1 and .psi..sub.2 to obtain a solution that converges
locally. This generalizes the Cauchy-Kovalevskaya extension
theorem.
Another form of this has slightly better convergence properties.
Starting from Equation (27) we write,
.mu..function..intg..times.
.psi..times..times..mu..function..psi..times..times..psi.
##EQU00088## noting that the initial conditions require .mu.(0)=1.
If D.sub.c.ident.(D.psi..sub.3).sup.-1D.sub.12, then we can write
this as the integro-differential operator, (1+.intg.D.sub.c).mu.=1,
where .intg. g.ident..intg..sub.0.sup.s g(.psi..sub.3)d.psi..sub.3.
One can now argue that,
.mu..infin..intg. ##EQU00089## which is equivalent to repeated
applications of (28). We can now use this theorem to find the field
reflected by the incident wavefront by assuming that
.PHI.(.psi..sub.1,.psi..sub.2)=-.alpha.n{tilde over
(F)}.sub.0(n{tilde over (x)}n)n, for x.di-elect
cons..differential.R defined by .psi..sub.3=0, and where F.sub.0(x)
is the incident (transmitted) wave.
It is possible to extend the Kovalevskaya extension to actually
recover the surface current paravector K for a perfect conductor
and hence the reflected wave. Suppose we are given the reflected
wave F.sub.r, we can recover the outgoing surface component of this
wave for .psi..sub.3.gtoreq.0 via,
.times..ident..psi..fwdarw. ##EQU00090## We actually take the Hardy
projection on to the boundary, which requires that the limit be
non-tangential, meaning that the limit sequence exists in a cone
smaller than a half plane about it's limit.
Our extended Kovalevskaya extension recovers the original function
from the boundary, namely .mu.P.sub.+F.sub.r=F.sub.r. We now
involve the surface projection operators, Q.sub.++Q.sub.-=1.
.mu..function..times..times..mu..times..times..times..times..times..mu..t-
imes..times..times..times..times..mu..times..times..times..times..times..i-
nfin..mu..times..times..times..times..times. ##EQU00091## We used
here the electromagnetic continuity condition at the conductor
boundary which requires Q.sub.-P.sub.+(F.sub.r+F.sub.0)=0, where
F.sub.0 is the incident wave.
This can be combined with Equation (28) to obtain the recursive
formulas,
.mu..function..psi..times..times..times..times..times..times..mu..functio-
n..intg..times..mu..function..psi..times..times..psi. ##EQU00092##
Convergence can be guaranteed for "small" .mu.Q.sub.+P.sub.+, which
is possible as .psi..sub.3.fwdarw.0, for sufficiently smooth
F.sub.0.
The main obstacle for use appears to be the proliferation of
symbolic terms for larger values of q.
Basis Functions and HOSVD
For a set of basis functions to approximate the scattered field, it
is convenient to look the separable solutions of the Helmholtz
equation, (D.sup.2+k.sup.2)F=0. There are 11 orthogonal coordinate
systems that are separable by the Helmholtz operator. Several of
those coordinate systems were described in Analytic Scattering
Models. We examine the spherical coordinate system and the
associated spherical harmonics that are solutions of the Helmholtz
equation, as one possible example for generating appropriate basis
functions.
The Laplacian in spherical coordinate systems can be written
as,
.times..psi..times..function..theta..times..function..theta..differential-
..differential..times..differential..psi..differential..differential..diff-
erential..theta..function..theta..times..differential..psi..differential..-
theta..function..theta..times..differential..psi..differential..PHI.
##EQU00093## If we set .psi.=R(kr).THETA.(.theta.).PHI.(.PHI.),
into the Helmholtz equation we obtain three separate ordinary
differential equations:
.times..times..times..times..function..times..function..theta..times..tim-
es..theta..times..function..theta..times..times..THETA..times..theta..func-
tion..theta..times..THETA..times..PHI..times..PHI..times..PHI.
##EQU00094## From this we can obtain the general solutions,
R(r)=B.sub.1h.sub.n.sup.1(kr)+B.sub.2h.sub.n.sup.2(kr), where
h.sub.n.sup.1(kr) is the spherical Hankel function of the first
kind, representing an inward traveling wave and h.sub.n.sup.2(kr)
is the spherical Hankel function of the second kind, representing
an outward traveling wave. These functions satisfy some recursion
relations and can be expressed in terms of Bessel functions of
fractional order. Additionally we have,
.times..function..times..times..function..times..times..times..times..tim-
es..times. ##EQU00095##
The polar angle dependency is solved by,
.THETA.(.theta.)=P.sub.n.sup.m(cos(.theta.)), where
P.sub.n.sup.m(x) is the Legendre polynomial of the first kind.
.PHI.(.PHI.)=exp(mj.PHI.).
The angular functions are typically collated into a single set of
functions known as spherical harmonics
Y.sub.n.sup.m(.theta.,.PHI.),
.theta..PHI..times..times..times..pi..function..times..function..theta..t-
imes..times..PHI. ##EQU00096## The spherical harmonics are
orthogonal satisfying,
.intg..times..pi..times..PHI..times..intg..pi.''.theta..PHI..times..theta-
..PHI..times..function..theta..times..times..theta..delta.'.times..times..-
delta.'.times. ##EQU00097## where .delta..sub.m'm is the
dirac-delta function. These functions are also complete, allowing
the representation of any function on the sphere, f(.theta.,.PHI.)
as a series of spherical harmonics. The spherical harmonics are
often written as a function of the unit normal vector
##EQU00098## as Y.sub.n.sup.m({circumflex over (x)}), which is a
polynomial function of the components of {circumflex over (x)},
called harmonic polynomials.
It is possible to decompose the spherical green's functions
h.sub.0.sup.1(kr) and h.sub.0.sup.2(kr) into spherical waves
when
.times.>.times.>.times. ##EQU00099## We can then write,
.function..infin..times..times. <.times.>.times. ##EQU00100##
where j.sub.l(z) is the spherical Bessel function of the first kind
and the real part of h.sup.1(z) for real z.
Note that, r= {square root over
(r.sub.1.sup.2+r.sub.0.sup.2-2r.sub.1r.sub.0 cos(.theta.))}, where
r.sub.0=.parallel.x.sub.0.parallel., r.sub.1=.parallel.x.parallel.,
and x.sub.0.sup.Tx=r.sub.1r.sub.0 cos(.theta.).
Suppose that r(r.sub.0,r.sub.1,cos(.theta.))= {square root over
(r.sub.1.sup.2+r.sub.0.sup.2-2r.sub.1r.sub.0 cos(.theta.))}, and
that we seek to find a decomposition similar to Equation (29) of
the form,
.function..function..theta..sigma..times..alpha..times..beta..times..gamm-
a..function..theta. ##EQU00101## Furthermore we seek such a
decomposition that gives us the best fit to r( ), with the minimum
number of terms. This problem is very similar to singular value
decomposition for matrices, but here the problem is multi-linear
containing multiple "eigen-vectors", .alpha..sub.l(r.sub.0),
.beta..sub.l(r.sub.1), and .gamma..sub.l(cos(.theta.)). This
problem can be solved using a device called multilinear singular
value decompositions, or higher order svd (HOSVD). This
decomposition reduces a multilinear form (tensor) into a sum of
mixed eigenvectors of the form,
.times..times..times..times..times..times..times..times..times..times..t-
imes..times..times..times. ##EQU00102## where the U.sub.ji.sup.(n)
are unitary matrices (entry (j,i)), and where
s.sub.i.sub.1.sub.i.sub.2.sub.. . . i.sub.N have all-orthogonal
subtensors and .parallel.s.sub.i.sub.1.sub.i.sub.2.sub.. . .
i.sub.N.parallel..sub.ij=q.sup.2 decreases as q increases from 1 to
I.sub.j. This means s.sub.i.sub.1.sub.i.sub.2.sub.. . . i.sub.N is
treated as a flat vector while holding one of the indices i.sub.j
fixed. The entries s.sub.i.sub.1.sub.i.sub.2.sub.. . . i.sub.N can
be sorted from largest to smallest norm, and the indices flattened
to make the decomposition look like,
.times..times..times..times..times..times..times..times..times..times.
##EQU00103## The indices for these tensors all represent the values
of continuous parameters, that have been densely sampled.
An example of the use of the HOSVD on the spherical Green's
function on the left hand side of Equation (29) is shown in FIG.
11. We allow r.sub.1 to range from 2.2 to 5 meters in increments of
a centimeter, and r.sub.0 ranges from 0.2 to 2.0 meters, while
cos(.theta.) is allowed to range from -1 to 1 in increments of
0.005. After applying the HOSVD, then flattening and sorting the
singular values, we plot the sorted singular values in FIG. 11.
Only 0.025 percent of the singular values are significant
contributors to the total value of the waveform. This means that a
massive simplification and compression is possible should one
tabulate the Green's function and decompose it into a separable
product of functions of r.sub.1, r.sub.0 and cos(.theta.).
Furthermore those functions of r.sub.1, r.sub.0 and cos(.theta.)
are smooth as can be seen by the locus plots in FIG. 12 and can be
approximated by differentiable splines. This greatly facilitates
the use of the integral formulas for the reproducing kernels in
Equation (20) and Equation (24). One uses a coordinate system in
this case whose origin is at the center of the scatterer and then
uses a kernel decomposition such as the HOSVD or the analytic
Equation (29) to evaluate the integral over the surface in closed
form. For the latter analytic case, the reproducing kernel itself
can be found by applying the D-kj operator to the left of each term
of the Green's function. Ideally however it would be preferable to
perform the Green's function decomposition using the coordinate
system appropriate for the reflector boundary type.
Harmonic functions generated by using HOSVD operations can be used
as basis functions in general. Indeed we can apply the Helmholtz
equation as a constraint, while matching Dirichlet style boundary
conditions for any given surface. Suppose we decomposed the
reproducing Kernel Equation (23) in a form similar to Equation
(30),
.times..times..times. ##EQU00104## where
x=xe.sub.1+ye.sub.2+ze.sub.3, is the receiver location and d.sub.l
are Clifford numbers. The field values are evaluated at the
receiver location x.sub.0, which is fixed for each receive antenna.
We can either use the reproducing kernel to produce the field at
the receiver, or we can use these functions to fit the boundary
conditions on a closed surface, using any convenient
parameterization of that surface. Typically we would place the
coordinate center at the center of the scattering object.
For the reproducing kernel case we have from Equation (20) and
Equation (32) for rectangular coordinates
.function..times..intg..differential..times..times..times..times.
##EQU00105## The scalar voltage could then be written as,
.intg..differential..times..times..times..times..times.
##EQU00106## where b is the vector effective length for the receive
antenna we are analyzing and F.sub.r(x) is found from Equation
(18). We can also use any coordinate parameterization for which the
scatterer represents a level surface, not just a rectangular
decomposition.
If 3-space is parameterized by .psi..sub.1,.psi..sub.2,.psi..sub.3,
and the closed surface of the scatterer is described by
.psi..sub.1=.rho., a constant we can use the HOSVD to decompose the
reflected wave
.function..psi..psi..psi..ident..times..psi..psi..psi..function..psi..psi-
..psi..apprxeq.
.times..alpha..psi..times..beta..psi..times..gamma..psi.
##EQU00107## The HOSVD operation is performed over a discretized
function of several variables, but we can use splines to
interpolate the function values .alpha..sub.l(.psi..sub.1),
.beta..sub.l(.psi..sub.2) and .gamma..sub.l(.psi..sub.3). We also
modify the algorithm to satisfy the Helmholtz equation,
(D.sup.2+k.sup.2)v(.psi..sub.1,.psi..sub.2,.psi..sub.3)=0 where
.psi..differential..differential..psi..times..psi..differential..differen-
tial..psi..times..psi..times..differential..differential..psi.
##EQU00108## The operator D.sup.2 will be a scalar operator that
can be discretized and treated as a multi-linear tensor contraction
on v(.psi..sub.1,.psi..sub.2,.psi..sub.3). This can be flattened to
simply look like a matrix. To obtain the approximation we attempt
to solve the Lagrangian problem,
.lamda..times..times..alpha..beta..gamma.
.times..intg..psi..times..psi..times..psi..times..function..rho..psi..psi-
..times.
.times..alpha..function..rho..times..beta..function..psi..times..-
gamma..function..psi..times..times..lamda..times..times..times.
.times..alpha..function..psi..times..beta..function..psi..times..gamma..f-
unction..psi. ##EQU00109## Here the integration
.intg..sub..psi.1.psi.2.psi.3 is approximated by a sum over a
region of points likely to represent a receiver/scatterer location.
Similar to the power method for eigenvalue problems, we can use
alternating directions optimization to obtain functions and
coefficients for this approximation problem.
Now implicitly v(.psi.) is dependent on the center of the
scattering body x.sub.r, which is one of the parameters we need to
estimate. After solving the Lagrange constrained minimization
problem in Equation (33) we obtain an approximating series solution
.SIGMA..sub.l
.sub.l.alpha..sub.l(.psi..sub.1).beta..sub.l(.psi..sub.2).gamma..sub.l(.p-
si..sub.3) also dependent on x.sub.r. Thus the voltage seen at
receiver position x.sub.0 is given by
V.sub.0(x.sub.r)=.SIGMA..sub.l
.sub.l.alpha..sub.l(.psi..sub.1.sup.0).beta..sub.l(.psi..sub.2.sup.0).gam-
ma..sub.l(.psi..sub.3.sup.0), where
(.psi..sub.1.sup.0,.psi..sub.2.sup.0,.psi..sub.3.sup.0) is the
parameterized coordinates of x.sub.0 relative to x.sub.r. If we
compute this over a dense set of possible reflector positions we
obtain another non-linear function that can be decomposed as a
function of x.sub.r=(x.sub.r,y.sub.r,z.sub.r). This function can be
decomposed again using HOSVD techniques to obtain,
.function..times.
.times..function..times..function..times..upsilon..function.
##EQU00110## Of course the voltage could theoretically be obtained
using any electromagnetic modeling tool, that can model the
scatterer and obtain the voltage at the receive antennas. The key
point here as that these computations can be done offline as part
of the calibration process. Only the final "tabulated" voltage as a
function of scatterer position needs to be retained for our
algorithm purposes.
Image Techniques
Similar to the electrostatic case, it is possible to place image
sources to satisfy electromagnetic boundary conditions. For a given
dipole moment vector b we can write the electromagnetic vector for
the simplest case of the Hertzian dipole as, F=E+icB, where,
.times..times..pi..times..times..times..times..times..times..times..pi..t-
imes..times..times..times..perp..function..times..times..times..function..-
times. ##EQU00111## where e.sub.r.ident.Dr=x/r and P
(e.sub.r)b=e.sub.r(e.sub.r b)=1/2(b+e.sub.rbe.sub.r) and
P.perp.(e.sub.r)b=b-e.sub.r(e.sub.r b)=1/2(b-e.sub.rbe.sub.r).
Suppose we have a Hertzian dipole 1 a distance h above an
electromagnetic screen. We intend to place an image dipole 2 at a
distance h below the screen, so that we can generate a reflected
field, -n{tilde over (F)}(n{tilde over (x)}n)n, whose superposition
with the incident field satisfies the boundary conditions.
This situation is illustrated in FIG. 13.
If the coordinate origin is on the electromagnetic screen, and
x.sub.1 is the center position of Dipole number 1 and n is the unit
normal to the conducting surface, then the coordinate center for
dipole 2, x.sub.2 is it's reflection through the conducting plane
P, ie x.sub.2=-nx.sub.1n. Note that the reflection operator is a
homomorphism, ie if R.sub.n(x).ident.n{tilde over (x)}n, then
R.sub.n(ab)=R.sub.n(a)R.sub.n(b), and R.sub.n(x)=x, when x.di-elect
cons.P. We can write the reflected field as,
.times..times..function..times..times..times..times..function..function.
##EQU00112## Suppose we decompose the E and B field into it's
tangential and normal components, E=E.sub.n+E.sub.t,
B=B.sub.n+B.sub.t, then we have, R.sub.n(E)=-E.sub.n+E.sub.t
R.sub.n(B)=-B.sub.n+B.sub.t. From this and the fact that the radii
are the same for both the reflected field and the original field,
we see that at boundary surface the total field is,
.times..times..times..times. ##EQU00113## This shows that the
correct boundary conditions are satisfied inasmuch as the total
tangential component of the electric field and the total normal
component of the magnetic field vanish.
Because of the homomorphism properties of the reflection operator
it is seen that the reflected field is the same as that due to an
image dipole with dipole moment b.sub.2=nbn at the reflected
position x.sub.2=-nx.sub.1n. The reflected electromagnetic field
becomes,
.times..times..pi..times..times..times..times..times..times..times..times-
..pi..times..times..times..times..perp..function..times..times..times..fun-
ction..times..times..times. ##EQU00114## ##EQU00114.2##
Imaging techniques are a powerful way to handle rectangular
boundaries. It's a good way to model electromagnetic propagation
near the ground or over water by treating the ground as an
electromagnetic conducting screen. Multiple plane boundaries can be
modeled by reflected each additional image across the new plane
introduced into the environment. In the next section we will show
how the concept can be extended to spherical boundaries using
conformal mappings in 4-space. Furthermore it is straightforward to
extend the technique to arbitrary spherical harmonics using the
reflection operator R.sub.n(F).
Conformal Mappings
A conformal map in a Clifford algebra takes the form of a Mobius
transform, .PHI.(x)=(ax+b)(cx+d).sup.-1. Where a, b, c and d are
members of the Clifford group .GAMMA..sub.3, which are products of
non-zero vectors in Cl.sub.0,3. These mappings preserve the surface
integrals in (16), namely
.intg..differential..times..function..times..times..function..intg..diffe-
rential..times..function..PHI..function..times..dagger..function..PHI..tim-
es..times..function..PHI..times..function..PHI..function.
##EQU00115## where
.function..PHI..ident..dagger. ##EQU00116## This implies,
.intg..times..function..times..function..times..times..intg..PHI..functio-
n..differential..times..function..PHI..function..times..dagger..function..-
PHI..times..times..function..PHI..times..function..PHI..function.
##EQU00117## Assuming h(y) is the reproducing kernel in (21), and
that g(y)) is left monogenic (ie satisfies the Helmoltz equation
(D+jk)g=0, we obtain,
.function..intg..PHI..function..differential..times..function..PHI..funct-
ion..PHI..function..times..dagger..function..PHI..times..times..function..-
PHI..times..function..PHI..function. ##EQU00118## where
y.sub.0=.PHI.(x.sub.0). This demonstrates that one can use
conformal maps to change the surfaces over which a solution can be
evaluated.
In particular the Mobius transform is known to map spheres to
spheres, where a "sphere" of infinite radius is a half plane. Thus
we can formulate a solution to the boundary value problem in the
half plane, using say the imaging solution of the prior section and
then use a conformal map to map the problem to a sphere or
vice-versa.
Another way we can use conformal maps is to use Vekua operator to
map a solution of the Laplace equation (k=0) to a solution of the
full Helmholtz equation. Conformal maps preserve the analytic
property of annihilating the Dirac operator. Thus if Dg(x)=0, then
DJ (.PHI.,x) g (.PHI.(x))=0. However we can use the Vekua operator
to map this to a solution that also satisfies the Helmholtz wave
equation D.sup.2u(x)+k.sup.2u(x)=0, via,
.function..function..PHI..times..function..PHI..function..intg..times..ti-
mes..times..times..function..times..times..times..function..PHI..function.-
.times. ##EQU00119## where J.sub.1( ) is the first order Bessel
function of the first kind. The Vekua transform is actually a
general way to associate solutions of the ordinary Laplace equation
with solutions of the Helmholtz equation.
Calibration
The calibration procedures for an array that is required to perform
geolocation consists of fitting a set of interpolating basis
functions to a set of observed voltages seen on the array, for
targets/reflectors at known positions in space. A lab or anechoic
chamber is typically used for initial calibration. Calibration
occurs not just for the free parameters of channel model, but also
can occur over different target types for a radar application, or
even dynamic targets that involve gestures, beating hearts, or
varying dielectric properties of materials.
Let v.sub.mk.sup.q(p.sub.q) be the voltage seen by the m'th sensor,
k'th frequency/channel and the q'th calibration collect, and where
p.sub.q is the "position" and shape parameter vector, for the q'th
calibration event. The calibration process attempts to fit a linear
combination of basis functions to the observed voltages. The fit
metric is usually least squares,
.mu..times..upsilon..function..times..alpha..times..times..times..functio-
n. ##EQU00120## The basis functions b.sub.lm(p.sub.q) are a
function of the position/shape parameter vector p.sub.q, antenna
index m and calibration index q. While it's generally not a good
ideal to calibrate in the presence of interference, we might also
presume that the received data and the basis functions have been
interference whitened over the antenna index. We of course expect
the basis functions to be different for different shapes in radar
applications whose signals are reflected off of scatterers of
different geometries.
The necessary conditions for the least squares calibration are
given by,
.times..times..times..function..times..times..alpha..times..times..times.-
.function..upsilon..function. ##EQU00121##
.times..times..times..times..function..times..times..times..function..tim-
es..alpha..times..times..times..function..times..upsilon..function.
##EQU00121.2## So for each antenna-frequency pair we have the
following L.times.L linear equation solution,
a.sub.mk=R.sub.b.sub.m.sub.b.sub.m.sup.-1R.sub.b.sub.m.sub.v.sub.mk
(37) where the l'th element of a.sub.mk is .alpha..sub.mk.sup.l the
l, n'th element of the matrix R.sub.b.sub.m.sub.b.sub.m is
.SIGMA..sub.qb*.sub.lm (p.sub.q)b.sub.nm (p.sub.q) and the l'th
element of R.sub.b.sub.m.sub.v.sub.mk is .SIGMA..sub.qb*.sub.lm
(p.sub.q) v.sub.mk.sup.q (p.sub.q).
It is interesting to note that certain coordinate systems allow the
interpolation of the wave equation over all space once the voltages
are defined over a closed surface. Thus is suffices to sample the
environment in the near-field over an appropriate closed surface.
In general however we are often required to grossly sub-sample the
array response. In this case we can often choose the few dominant
terms in our basis function, either found using HOSVD as in Basis
Functions and HOSVD, or using something similar to Generalized
Spherical harmonics from the same Section. The higher order
harmonics have decreasing contributions to the overall wave
function, especially in the far field.
Obtaining The Basis Functions
From Channel Modeling as described herein, a variety of techniques
for generating analytic formulas for the wavefunction are shown,
even for complex scatterers. Here we give a brief overview of the
procedure.
1. Generate a model for the transmitted electromagnetic field.
2. Predict the field being scattered (for radar applications) off
of the desired target type.
3. Use the prediction to generate basis functions for the received
voltage at each antenna receiver and frequency,
4. Collect data in the near field for multiple scatterer/receiver
positions,
5. Use least squares to estimate the complex combining weights for
the basis functions.
The model for the transmitted electromagnetic field is at least
partially determined by numerical modeling of the antenna array. We
have seen a formula for the Hertzian dipole given in Equation 35.
This formula actually corresponds to the first order Spherical
Harmonics. In general we might write the electromagnetic vector
potential as,
A.sub.u(x)=h.sub.n(u).sup.2(kr)Y.sub.n(u).sup.m(u)({circumflex over
(x)})b, where the u index has been flattened to include all m and n
components of the spherical Harmonics, and b is a dipole moment
vector of the transmitting antenna. We can also use the
h.sub.n(u).sup.1 (kr) for incoming waves. Similar considerations
can be made for other coordinate systems.
From this we obtain the electric field from Equation (22) as,
.function..times..times..PHI..function..function. ##EQU00122##
.PHI..function..times..function. ##EQU00122.2##
.function..times..function. ##EQU00122.3## We then use one of the
several scattering modeling techniques described in Analytic
Scattering Models to obtain an induced scattered field
.sub.u(x.sub.m;p.sub.q) at receiver position x.sub.m and scatterer
parameter vector p.sub.q. The scattered field itself may require a
series solution for each transmitted E.sub.u(x) and we may see a
multi-parameter basis function of the form,
.sub.pu(x.sub.m;p.sub.q), whose indices (p,u) can either be index
flattened to .sub.l(x.sub.m; p.sub.q), or further sub-selected on
the basis of the largest singular values using HOSVD techniques
described in Basis Functions and HOSVD.
Finally the actual voltage and the scalar basis function is
obtained via the vector effective length for antenna m,
b.sub.lm(p.sub.q)=b.sub.m .sub.l(x.sub.m;p.sub.q). Once the complex
gains for the model are determined from (37), we have a series
solution model for the received voltage that can be used in the
parameter estimation and tracking models of Geolocation and
Tracking.
Calibration Model Modifications
There are a number of modifications to the basis function model we
have proposed so far that may come into play. The first occurs when
we only have a sparse number of calibration points to work with.
This can happen when calibration is very expensive or can not be
automated.
We can relax the requirement of learning a new set of gains for
every frequency index k, by creating a structured version of our
calibration gains, namely:
.alpha..sub.mk.sup.l=.beta..sub.m.sup.l.gamma..sub.k. This asserts
that the basis functions via it's dependency on the wavenumber k
suffices to interpolate over frequency, outside some frequency
dependent gains on the transmit to receive antenna path. The
structured gains can be found by holding the .beta..sub.m.sup.l
terms constant and minimizing over .gamma..sub.k, followed by
holding the .gamma..sub.k constant and optimizing over
.beta..sub.m.sup.l. This is repeated until convergence.
Another modification occurs if we use the pre-computed grid points
suggested in (25). This would use more of a discrete point model
for the reflector surface points. Those values, rather than the
basis function coefficients become the output of the calibration
process, though modeling the transmit antenna response may well use
interpolating basis functions as described earlier.
Finally we may find that multiple nearby receive antennas mutually
couple. If we model these antennas as re-radiating scatterer, than
each antenna has a fixed contribution to another nearby antenna
that is a complex multiply of it's receive voltage. This can be
modeled as an additional mixing matrix added to the final received
voltage, ie,
.upsilon..function..apprxeq.'.times..times..times.'.times..times..alpha.'-
.times..times..times..times.'.function. ##EQU00123## The mixing
matrix can also be discovered as part of the calibration process.
Again we use iterative alternating directions least squares
optimization for each linear component as a technique for
determining these coefficients from a sufficiently dense
calibration table v.sub.mk.sup.q(p.sub.q).
In one or more of the embodiments of the present invention there is
provided a method of parameter estimation in a multi-channel signal
environment system where a plurality of receive antennas receive
signals from one or more transmitters that transmit a signal or
wave that is reflected from one or more targets, or receive
antennas that receive directly from the transmitters, whose
received signals are processed over multiple frequencies or
channels by a digital receiver. One or more processors may be
provided with processing comprises the steps of: (a) calibrating
before the operation of the digital receiver to determine the free
parameters of a mathematical model of a channel either as the
channel model parameters or in the form of tabulated data; (b)
calibrating during the operation of the digital receiver to
determine the channel model; (c) comparing received antennas array
voltages to an analytic or table driven channel model from a
calibrated template without only relying on information lossy
intermediate steps such as delay time of arrival or angle of
arrival measurements; (d) creating a statistical likelihood
function modeling receiver noise to determine model channel
parameters or prior channel uncertainty; (e) saving target
reflector/emitter parameters to be reused for dynamic tracking; and
(f) using Bayesian particle filtering or Maximum Likelihood Methods
to update mixture models for the unknown parameters.
The method may use Bayesian detection or other Statistical Signal
Processing Techniques for the estimation of channel parameters such
as location parameters, shape parameters and reflector
electromagnetic parameters. Various aspects of the invention could
be defined wherein:
the transmitted signal is a frequency stepped radar;
a statistical likelihood function is used to determine target type
or target position;
a channel template for each target type is used to further classify
the target type or estimate the target position;
static direct path clutter is removed via a cancellation
algorithm;
the parameter estimation blindly calculates an unknown gain and
phase constant .beta.q over each channel observation q, to absorb
bulk phase and gain changes between observations;
the statistical likelihood function takes the form of a Cauchy
quotient and is therefore a member of the numerical range of a
parameterized matrix as exhibited in Equation (14).
the target classification is determined from the likelihood
function using the posterior probability of detection via Bayes
theorem;
the target classification is determined from the likelihood
function after blindly estimating any position dependent parameters
during the classification process itself;
the target classification is determined from the likelihood
function after blindly estimating phase terms such as time delays
of arrival;
the target classification is determined from the likelihood
function after blindly estimating phase delays using fast Fourier
transform (FFT) processing;
the target classification or position parameters are determined
from the likelihood function that assumed the scale parameters
.beta.q=.beta., are independent of the collection index and can
therefore be solved using the generalized eigenvalue problem in
Equation (8);
the DSP processing hardware is configured to enable parallel
likelihood function calculations wherein the same instruction
across all processing resources;
the DSP hardware has built in support for parallel reduction for
primitive associative operations like addition, maximization,
minimization or multiplication to enable the parallel computation
of the likelihood function;
multiple transmit and receive RF chains are packed into an radio
frequency (RF) integrated circuit (IC) with clocks controlled by a
phased locked loop;
the signals from each antenna are down converted and digitized by a
bank of analog to digital converters, where they are further
decimated to obtain a sample of the channel for each antenna and
frequency in use;
the wave function can be further interpolated using basis functions
that are reflected through objects of predetermined shape or
presumed electrical properties;
interpolating uses a likelihood function to compute a fixed set of
channel vectors .alpha., modulating the basis functions, that can
be used to interpolate the channel at an arbitrary target
position;
the method further includes blindly calculating an unknown gain and
phase constant .beta.q over each channel observation q, to absorb
bulk phase and gain changes between observations;
the method further includes computing separate channel vectors
.alpha..sub.c for each target type c;
the method further includes jointly computing the phase constants
.beta.q and the channel vectors .alpha..sub.c by alternating
directions optimization; ie iteratively estimating each parameter
while holding the other fixed;
the method further includes assuming the inner product of the basis
functions is independent of position resulting in an eigenvalue
problem (10) that is used to learn the channel vectors during
calibration;
the method further includes assuming the scale parameters
.beta.q=.beta., are independent of the collection index, and can
therefore be solved using the generalized eigenvalue problem in
Equation (8);
the method further calibrates a blind classifier by blindly
estimating any position dependent parameters during the
classification process itself, without necessarily using known
target or array transceiver positions;
the method further calibratesg a blind classifier by after blindly
estimating phase terms such as time delays during the
classification process itself without necessarily using known
target or array transceiver positions;
the method further calibrates a blind classifier by blindly
estimating phase delays using FFT processing, during the
classification process itself, without necessarily using known
target or array transceiver positions;
the DSP processing hardware is configured to enable parallel
likelihood function calculations, yet can share the same
instruction across all processing resources;
the DSP hardware has built in support for parallel reduction for
primitive associative operations like addition, maximization,
minimization or multiplication to enable the parallel computation
of the likelihood function;
multiple transmit and receive RF chains are packed into an RF IC,
with clocks controlled by a phased locked loop;
the signals from each antenna are down-converted and digitized by a
bank of analog to digital converters wherein the signals are
further decimated to obtain a sample of the channel for each
antenna and frequency in use; OR
the likelihood function is itself used to further refine the
presumed positions of the transmitter and receiver to improve the
calibration channel model.
the foregoing and as mentioned above, it is observed that numerous
variations and modifications may be effected without departing from
the spirit and scope of the novel concept of the invention. It is
to be understood that no limitation with respect to the embodiments
illustrated herein is intended or should be inferred. It is
intended to cover, by the appended claims, all such modifications
within the scope of the appended claims.
* * * * *