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.

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 Estaimation 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,

$\begin{array}{cc}{x}_k\left(t\right)=\sum _{q=1}^Q{\mathscr{H}}_{\mathrm{kq}}\left({\alpha }_{\mathrm{kq}},p_q,s_{\mathrm{kq}}\left(t\right)\right)+i_k\left(t\right),& \left(1\right)\end{array}$

where x_k(t) is the M_k×1 received waveform vector at sensor k and time t, i_k(t) is an additive interference plus noise process for sensor k, _kq (α_kq,p,s_kq(t)) is a vector channel process that modulates the transmitted signal s_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, α_kq and a global parameter vector p_q, which is of primary interest to our receiver processing algorithms. The p_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,

$P\left(x\text{|}p\right)=\prod _{k=1}^K{\int }_{{\alpha }_k}f_k\left(x_k-\sum _{q=1}^Q{\mathscr{H}}_{\mathrm{kq}}\left({\alpha }_{\mathrm{kq}},p_q,s_{\mathrm{kq}}\left(t\right)\right)\right)q_k\left({\alpha }_{\mathrm{kq}}\right)d{\alpha }_{\mathrm{kq}},$

where P(x|p) is the probability of observing the total view column vector x=[x₁^T,x₂^T, . . . ,x_K^T]^T, given p=[p₁^Tp₂^T, . . . p_Q^T]^T, f_k (x) is the probability density function for the interference process i_k, and g_k (α_kq) is the prior distribution of the noisy, free parameters of the channel process _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,

$\begin{array}{cc}{x}_k\left(t\right)=\sum _{q=1}^Q{\mathscr{H}}_{\mathrm{kq}}\left({\alpha }_{\mathrm{kq}},p_q\right)s_{\mathrm{kq}}\left(t\right)+i_k\left(t\right).& \left(2\right)\end{array}$

The notation can be further simplified by writing this in the time sampled digital domain as,

x _k(n)=H_ks_k(n)+i_k(n),

H_k≡[_k1(α_k1,p₁),_k2(α_k2,p₂) . . . ,_kQ(α_kQ,p_Q)]

where n is the time sample index, H_k is the M_k×Q channel matrix as a function of the nuisance parameter vector α_k≡[α_k1^T, α_k2^T . . . α_kQ^T] and parameter of interest vector p≡[p₁^T,p₂^T, . . . p_Q^T]^T, for all Q emitters/reflectors. Also the emitter signal waveforms are packed into the Q×1 signal vector s_k(n)≡[s_k1(n),s_k2(n) . . . ,s_kQ(n)]^T.

Colored Gaussian Interference

Some processing simplifications can be made by subsuming all unmodeled interfering waveforms into the single interference vector i_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 _i(i_k|R_ii)=|πR_ikik|^-1 exp(−i^HR_ikik^-1i).

The zero mean Gaussian normal distribution is referred to by, N (0, R_ikik), and one can also define the complex normal distribution function as,

N(i,R)≡|πR|^-1 exp(−i^HR^-1i).

This allows the likelihood function to be written as,

$p\left(x_kH_k,s_k\right)=\prod _{n=1}^N{\pi R_{i_ki_k}}^{-1}\exp \left(-\mathrm{Tr}\left(R_{i_ki_k}^{-1}o\left(x_x\left(n\right)-H_ks_k\left(n\right)\right)\right)\right),$

where Tr(⋅) is the matrix trace operator and o(a)≡aa^H is the outer product operator.

We now perform some manipulations to expose the sufficient statistics for the probability function as follows:

$\begin{array}{cc}\begin{array}{c}p\left(x_kH_k,s_k\right)={\pi R_{i_ki_k}}^{-N}\exp \left(-\mathrm{Tr}\left(R_{i_ki_k}^{-1}\sum _{n=1}^No\left(x_k\left(n\right)-H_ks_k\left(n\right)\right)\right)\right)\\ ={\pi R_{i_ki_k}}^{-N}\mathrm{etr}(-{\mathrm{NR}}_{i_ki_k}^{-1}({\hat{R}}_{x_kx_k}-H_k{\hat{R}}_{x_ks_k}^H-\dots \\ {\hat{R}}_{x_ks_k}H_k^H+H_k{\hat{R}}_{s_ks_k}H_k^H)),\end{array}\hspace{1em}& \left(3\right)\end{array}$

where the time averaged cross and auto correlation matrices are defined by

${\hat{R}}_{\mathrm{AB}}\equiv \frac{1}{N}\sum _{n=1}^NA\left(n\right)B^H\left(n\right),$

and where etr (A)≡exp (Tr(A)). A matrix version of completing the square provides a way to identify sufficient statistics for this problem. We write:

p(x_k|H_k,s_k)=|πR_ikik|^-N etr(−NR_ikik^-1({circumflex over (R)}_xkxk−{circumflex over (R)}_xksk{circumflex over (R)}_sksk^-1{circumflex over (R)}_xksk^H+ . . . (H_k−{circumflex over (R)}_xksk{circumflex over (R)}_sksk^-1){circumflex over (R)}_sksk(H_k−{circumflex over (R)}_xksk{circumflex over (R)}_sksk^-1)^H)).

Note that Ĥ_k≡{circumflex over (R)}_xksk{circumflex over (R)}_sksk^-1 is a least squares approximation of the channel matrix, and that {circumflex over (R)}_ikik≡{circumflex over (R)}_xkxk−{circumflex over (R)}_xksk{circumflex over (R)}_sksk^-1{circumflex over (R)}_xksk^H is an estimate of the interference covariance matrix. This allows one to write,

p(x_k|H_k,s_k)=η etr(−NR_ikik^-1((H_k−Ĥ_k){circumflex over (R)}_sksk(H_k−Ĥ_k)^H))   (4)

η≡|πR_ikik|^-N etr(−NR_ikik^-1{circumflex over (R)}_ikik).

This makes Ĥ_k≡{circumflex over (R)}_xssk{circumflex over (R)}_sksk^-1 a sufficient statistic given R_ikik and s_k.

We can write this in the whitened form,

p(x_k|H_k,s_k)=η etr(−N(Ĥ_sk−{circumflex over ({tilde over (H)})}_sk)^H({tilde over (H)}_sk−{circumflex over ({tilde over (H)})}_sk))   (5)

where {tilde over (H)}_sk≡R_ik^-HH_k{circumflex over (R)}_sk^H is the interference whitened and signal power normalized channel matrix for view k, and where {dot over ({tilde over (H)})}_sk≡R_ik^-H{circumflex over (R)}_xksk{circumflex over (R)}_sk^-1, is an estimate of the interference whitened and signal power normalized channel matrix. The matrix R_ik is the upper triangular Cholesky factor of R_ikik with R_ik^-H being the Hermitian transpose of the inverse of that matrix. Similarly {circumflex over (R)}_sk^H is the Hermitian transpose of the Cholesky factor of {circumflex over (R)}_sksk.

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_kq(n) being a sufficiently wide band signal to make {circumflex over (R)}_sksk a full rank matrix. We must also cover the very important case of {circumflex over (R)}_sksk 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_k is actually a rank-1 vector modeled as the sum of it's multipath components,

$\begin{array}{cc}{H}_k\equiv h_k=\sum _{q=1}^Qh_{\mathrm{kq}}\left(p\right){\alpha }_{\mathrm{kq}}.& \left(6\right)\end{array}$

Here we assume that the position parameters p is shared by all reflected multipath components.

Narrow Band Channels

The important case wherein s_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_k and q is associated with a path delay τ_q, we can approximate the signal of interest as,

s_kq(n)≈s_k(n)e^-2πjfkτ_q,

where s_k (n) is a narrow band signal centered at f_k.

{circumflex over (R)}_xksk={circumflex over (r)}_xkskv_k^H,

where

${\hat{r}}_{x_ks_k}=\frac{1}{N}\sum _{n=1}^Nx_k\left(n\right)s_k^{*}\left(n\right),\mathrm{and}v_k\equiv {\left[e^{-2\pi {\mathrm{jf}}_k{\tau }_1},e^{-2\pi {\mathrm{jf}}_k{\tau }_2}\dots e^{-2\pi {\mathrm{jf}}_k{\tau }_Q}\right]}^T.$

Similarly we have,

{circumflex over (R)}_sksk≡{circumflex over (r)}_skskv_kv_k^H.

This allows Equation (3) to be written as,

$\begin{array}{c}p\left(x_kH_k,s_k\right)={\pi R_{i_ki_k}}^{-N}\mathrm{etr}(-{\mathrm{NR}}_{i_ki_k}^{-1}({\hat{R}}_{x_kx_k}-H_kv_k{\hat{r}}_{x_ks_k}^H-\dots \\ {\hat{r}}_{x_ks_k}v_k^HH_k^H+{\hat{r}}_{s_ks_k}H_kv_kv_k^HH_k^H)),\\ ={\pi R_{i_ki_k}}^{-N}\mathrm{etr}(-{\mathrm{NR}}_{i_ki_k}^{-1}({\hat{R}}_{i_ki_k}-\dots \\ \left({\hat{h}}_k-h_k\right){{\hat{r}}_{s_ks_k}\left({\hat{h}}_k-h_k\right)}^H)),\end{array}$ $\mathrm{where}{\hat{h}}_k\equiv {\hat{r}}_{x_ks_k}/{\hat{r}}_{s_ks_k},h_k=H_kv_k,\mathrm{and}$ ${\hat{R}}_{i_ki_k}\equiv {\hat{R}}_{x_kx_k}-{\hat{r}}_{x_ks_k}{\hat{r}}_{x_ks_k}^H/{\hat{r}}_{s_ks_k}.$

Similar to Equation (5), one can write the whitened version of the likelihood function for the narrow band channel as,

p(x_k|H_k,s_k)=η exp(−N{circumflex over (r)}_sksk∥{tilde over (h)}_k−{tilde over (ĥ)}_k∥²)   (7)

η≡|πR_ikik|^-N etr(−NR_ikik^-1{circumflex over (R)}_ikik),

where {tilde over (h)}_k≡R_ik^-Hh_k and {tilde over (ĥ)}_k≡R_ik^-Hĥ_k.

Aperture Noise Model

For proper estimation of channel parameters it is critical to have an accurate noise model of the aperture vector _kq (α_kq,p_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=_kq (α_kq, p_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_t,

h=h _t +e,

where e is colored Gaussian complex noise. In general e may be dependent on the unknown parameters α_kq, p_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 _t⊙(1+e),

where ⊙ 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_t⊙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_k×1 vector,

{tilde over (h)}_kα_k,

where {tilde over (h)}_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)} _k ={tilde over (h)} _kα_k+e_k

α_k=1+ε_k.

Here, {tilde over (x)}_k is the whitened received data, in the form of an estimated channel response, playing the role of {tilde over (ĥ)}_k in Equation (7).

Assuming a complex Gaussian noise model we have,

$\begin{array}{cc}P\left({\tilde{x}}_k,{\alpha }_k{\tilde{h}}_k\right)=\frac{1}{{\pi }^{M_k}{\sigma }^{2M_k}{\beta }_k^2}\exp \left(-\frac{1}{{\sigma }^2}{{\tilde{x}}_k-{\tilde{h}}_k{\alpha }_k}^2-\frac{1}{{\beta }_k^2}{{\alpha }_k-1}^2\right).& \left(8\right)\end{array}$

Here σ² is the variance of the additive white Gaussian noise (post whitening) and β_k² is the variance that controls the deviation of the per-frequency channel gain from unity. By varying β_k² we can express our confidence in the fidelity of the channel model. A large β_k² 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 α_k to obtain,

P({tilde over (x)}_k|{tilde over (h)}_k)=det(πR_xkxk)^-1 exp(−({tilde over (x)}_k−{tilde over (h)}_k)^HR_xkxk^-1({tilde over (x)}_k−{tilde over (h)}_k))   (9)

where

R_xkxk≡σ²I+{tilde over (h)}_k{tilde over (h)}_k^Hβ_k².

The posterior expectation for α_k and it's maximum likelihood estimate is given by,

${\hat{\alpha }}_k=\frac{{\beta }_k^2{\tilde{h}}_k^H{\tilde{x}}_k+{\sigma }^2}{{\beta }_k^2{\tilde{h}}_k^H{\tilde{h}}_k+{\sigma }^2}.$

It can be shown that to evaluate the argument in the exponential of Equation (9), it suffices to substitute {circumflex over (α)}_k into α_k in Equation (8). Also the determinant is easily computed as,

det(πR_xkxk)=π^Mk(σ²+β_k²∥{tilde over (h)}_k∥²).

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)} _mk ={tilde over (h)} _mkα_mk+e_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)}_k)=det(πR_xkxk)^-1 exp(−({tilde over (x)}_k−{tilde over (h)}_k)^HR_xkxk¹({tilde over (x)}_k−{tilde over (h)}_k)),

where

R_xkxk≡σ²I+(|{tilde over (h)}_k|²)(β_k²),

where (|{tilde over (h)}_k|²) is an M×M diagonal matrix, whose m′th diagonal element is |{tilde over (h)}_mk|² and (β_k²) is a diagonal matrix whose m′th diagonal element is β_mk², the variance of the gain multiplier α_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 _k =H _k g _k +e _k   (10)

H _k =[h _k1 ,h _k2 , . . . h _kQ]

g_k≡[α_k1,α_k2, . . . α_kQ]^T

g_k˜N(1,R_gkgk),

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_k˜N(H_k1, R_ekek+H_kR_gkgkH_k^H), so that we can write the likelihood function as,

p(x_k|H_k)=det(πR_xkxk)^-1 exp(−(x_k−H_k1)^HR_xkxk^-1(x_k−H_k1))

R _{x k x k} =R _{e k e k} +H _k R _{g k g k} H _k ^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,

$x_k=a_k+u_k,a_k=\sum _{q=1}^Qh_{\mathrm{kq}}{\alpha }_{\mathrm{kq}}=H_kg_k,$

where the target channel α_k is seen as a mixture of multiple related component channels, (e.g. a human body with torso and limbs). The u_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_k˜N(g_k,R_gg)

u_k˜N(ū_k,R_uu).

The probability for the signal on and signal off are,

p(x_k|ON)=N(x_k−H_kg_k−ū_k,H_kR_ggH_k^H+R_uu)

p(x_k|OFF)=N(x_k−ū_k,R_uu),

from this one can compute the posterior probability of the target present (ON) and target not present (OFF),

$p\left(\mathrm{ON}x_k\right)=\frac{p\left(x_k\mathrm{ON}\right)p_1}{p\left(x_k\mathrm{ON}\right)p_1+p\left(x_k\mathrm{OFF}\right)p_0},$

where p₁ is the prior probability of the target being present and p₀=1−p₁.

An additional modification to our model is to allow the removal of bulk phase terms for each channel component. Thus we have,

$x_k=\sum _{q=1}^Qh_{\mathrm{kq}}{\alpha }_{\mathrm{kq}}c_q+u_k,$

where c_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_q, simply by rewriting,

$x_k-\sum _{q\ne p}h_{\mathrm{kq}}{\alpha }_{\mathrm{kq}}c_q-u_k=h_{\mathrm{kp}}c_p$ ${\hat{c}}_p=\mathrm{sgn}\left(h_{\mathrm{kp}}^H\left(x_k-\sum _{q\ne p}h_{\mathrm{kq}}{\alpha }_{\mathrm{kq}}c_q-u_k\right)\right),\mathrm{sgn}\left(u\right)\equiv \frac{u}{u}.$

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,

$h_p\approx x_k-\sum _{q\ne p}h_{\mathrm{kq}}{\alpha }_{\mathrm{kq}}c_q-u_k.$

We either pick h_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_kq (p_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_q may contain the uncertain scatterer/emitter positions. We may assume or estimate priors for these parameters, p_q˜f(p_q) and u˜g(u). If these prior distributions are Gaussian, then we have, p_q˜N (p_0q;R_p) and u˜N (u₀,R_u) . Defining the real multivariate normal distribution by (α,R)≡det(2πR)^-1/2exp(−½α^TR^-1α), we can write the likelihood function as,

$p\left(x_k,p_0,p,u\right)=p\left(x_kH_k\right)\left(u-u_0,R_u\right)\prod _{q=1}^Q\left(\left(p_q-p_{0q},R_p\right)p\left(p_{0q}\right)\right).$

Using Bayes theorem we can write out the posterior distribution as,

$\begin{array}{cc}p\left(p_0x_k\right)=\frac{{\int }_{p_1}{\int }_{p_2}\dots {\int }_{p_Q}{\int }_up\left(x_k,p_0,p,u\right)}{{\int }_{p_1p_{01}}{\int }_{p_2p_{02}}\dots {\int }_{p_Qp_{0Q}}{\int }_up\left(x_k,p_0,p,u\right)}.& \left(11\right)\end{array}$

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 _km =h _km e ^{-2πjf k τ m} +i _km,

where τ_m is the signal delay seen at antenna m, and where f_k is the signal frequency for index k. We perform an inverse discrete Fourier transform, ^-1 and obtain,

${}^{-1}\left(x_{\mathrm{km}}\right)=\sum _{k=0}^{K-1}e^{2\pi {\mathrm{jf}}_kt}x_{\mathrm{km}}=\sum _{k=0}^{K-1}h_{\mathrm{km}}e^{2\pi {\mathrm{jf}}_k\left(t-{\tau }_m\right)}+\sum _{k=0}^{K-1}e^{2\pi {\mathrm{jf}}_kt}i_{\mathrm{km}}.$

If its assumed that the frequency hopping is linear f_k=f₀+kΔ, and that h_km≈h_m is dominated by a frequency independent mode, we then have,

$\begin{array}{c}{}^{-1}\left(x_{\mathrm{km}}\right)=h_me^{2\pi {\mathrm{jf}}_0\left(t-{\tau }_m\right)}\frac{e^{2\pi \mathrm{jK}\Delta \left(t-{\tau }_m\right)}-1}{e^{2\pi j\Delta \left(t-{\tau }_m\right)}-1}+{\stackrel{⋓}{i}}_{\mathrm{km}}\\ =h_me^{2\pi j\left(\left(f_0+\frac{K-1}{2}\Delta \right)\left(t-{\tau }_m\right)\right)}\frac{\sin \left(\frac{K\Delta }{2}\left(t-{\tau }_m\right)\right)}{\sin \left(\frac{\Delta }{2}\left(t-{\tau }_m\right)\right)}+{\stackrel{⋓}{i}}_{\mathrm{km}}.\end{array}$

One can plot what the simple delay channel looks like in the time domain. Here assuming K=128, Δ=20 MHz and τ_m=5 ns, we plot the impulsive response in FIG. 5.

As the number of frequencies K→∞, 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,

$\begin{array}{cc}p\left(p\text{|}x\right)=\frac{p\left(x\text{|}p\right)p\left(p\right)}{\int p\left(x\text{|}p\right)p\left(p\right)\mathrm{dx}},& \left(12\right)\end{array}$

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),

$p\left(x\text{|}p\right)=\prod _{k=1}^K{\det \left({\mathrm{\pi R}}_{x_kx_k}\left(p\right)\right)}^{-1}\exp \left(-{\left(x_k-h_k\left(p\right)\right)}^HR_{x_kx_k}^{-1}\left(p\right)\left(x_k-h_k\left(p\right)\right)\right),\mathrm{where},R_{x_kx_k}\left(p\right)\equiv {\sigma }^2I+h_k\left(p\right)h_k^H\left(p\right){\beta }_k^2.$

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πR) ^-1/2 exp(−½p^TR^-1p).

One can propose to model the prior and posterior probabilities as a Gaussian mixture. The prior is assumed to adhere to,

$p_0\left(p\right)=\sum _{j=1}^J{\eta }_j\left(p-m_j;R_j\right),$

where Σ_jη_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₀(p).

Gaussian Mixture Filter Update

1. Sample N random points from the Gaussian mixture prior p₀(p), p_n.

2. Compute the conditional likelihood function over all the random points,

λ_n=p(x|p_n).

3. Normalize the likelihoods to sum to one,

${\lambda }_n\to {\lambda }_n\text{/}\sum _n{\lambda }_n$

4. Resample a subset of Q<N points from the original N, where the probability of choosing p_n equals λ_n. Call this set of renumbered points p_v(q) for subindex q and the associated original index, n=v(q).

5. Use the EM algorithm, and the points p_v(q) to learn a new set of n′_j, m′_j and R′_j to fit the set of sampled p_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,

$p_0\left(p\right)=\sum _{j=1}^J{\eta }_jf_j\left(p\right),$

we can generate N pseudo-random samples from any such distribution for which we can sample directly from f_j(p). For each n we randomly select what category j, p_n will belong to, based on the category probabilities η_j. Selecting the category can be done in approximately o(log(J)) comparisons by using a binary search into the cumulative sum array of η_j, namely ç_j≡Σ_q=1^jη_j. When a uniform random variate u is chosen, we need merely find j such that, ç_j-1≤u<ç_j, using binary search, and initializing ç₀=0.

Once the category is chosen for a given n, we generate a random sample using the pdf f_j(p). For example when f_j(p)=(p−m_j;R_j), we can sample from any given multivariate Gaussian distribution by setting

p _m =m _j C _j ^T e _n,

where e_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_j matrix is the upper Cholesky factor of the positive definite covariance matrix R_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

${\lambda }_n=\frac{p\left(x\text{|}p_n\right)}{{\sum }_np\left(x\text{|}p_n\right)}.$

Again we can compute a cumulative sum of the likelihoods v_n≡Σ_q=1ⁿλ_n, followed by a binary search to find a n where a given, newly generated random uniform variable u_q satisfies, v_n−1≤u<v_n. This is our n=v(q).

Let us compute the probability distribution for p_v(q). We have,

$\begin{array}{c}p\left({p_v}_{\left(q\right)}=p\text{|}x\right)=\sum _{n=1}^Np\left(p\text{|}n=v\left(q\right),x,p_n=p\right)p\left(n=v\left(q\right)\text{|}x,p_n=p\right)\\ =\sum _np_0\left(p\right){\lambda }_n\text{|}\left(x,p_n=p\right)\\ =\sum _{n=1}^N\frac{p\left(x\text{|}p\right)p_0\left(p\right)}{{\sum }_{q\ne n}p\left(x\text{|}p_q\right)+p\left(x\text{|}p\right)}\\ \approx \frac{{\sum }_{n=1}^Np\left(x\text{|}p\right)p_0\left(p\right)}{N{\int }_pp\left(x\text{|}p\right)p_0\left(p\right)\mathrm{dp}}\\ =\frac{p\left(x\text{|}p\right)p_0\left(p\right)}{{\int }_pp\left(x\text{|}p\right)p_0\left(p\right)\mathrm{dp}}\\ =p\left(p\text{|}x\right).\end{array}$

We have Σ_q≠np(x|p_q)+p(x|p)≈N ∫_pp(x|p)p₀(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_v(q) asymptotically approaches our desired posterior distribution. In what follows we label the resampled points as p_v(q)→p_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,

$p_1\left(p\text{|}m,R\right)=\sum _{j=1}^J{\pi }_j\left(p-m_j;R_j\right).$

The initial values for π_j, m_j and R_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₁(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_j and R_j, to find the next iterate. This can be summarized as,

Expectation Step

Compute the posterior likelihood of a sampled point p_q being in category j,

$\begin{array}{cc}{\gamma }_{\mathrm{qj}}=\frac{\pi j\left(p_q-m_j,R_j\right)}{{\sum }_{k=1}^J{\pi }_k\left(p_q-m_k,R_k\right)}N_j\equiv \sum _{q=1}^Q{\gamma }_{\mathrm{qj}}.& \left(13\right)\end{array}$

The parameter N_j can be thought of as the expected number of points associated to category j. Note Σ_jN_j=Q.

Maximization Step

The updated mixture parameters are given by,

$m_k^{\mathrm{new}}=\frac{1}{N_k}\sum _{q=1}^Q{\gamma }_{\mathrm{qk}}p_q$ $R_k^{\mathrm{new}}=\frac{1}{N_k}\sum _{q=1}^Q{\gamma }_{\mathrm{qk}}\left(p_q-m_k^{\mathrm{new}}\right){\left(p_q-m_k^{\mathrm{new}}\right)}^T$ ${\pi }_k^{\mathrm{new}}=\frac{N_k}{Q}.$

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,

$L=\sum _{q=1}^Q\log \left(\sum _{j=1}^J{\pi }_j\left(p_q-m_j;R_j\right)\right).$

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_n drawn from the prior distribution, modify (13) to read,

${\gamma }_{\mathrm{nj}}=\frac{{\pi }_j\left(p_n-m_j,R_j\right)p\left(x\text{|}p_n\right)}{{\sum }_{k=1}^J{\pi }_k\left(p_n-m_k,R_k\right)p\left(x\text{|}p_n\right)}.$

This causes a given p_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 ρ(p|x) is the posterior distribution of position parameter p given the observed data x, after many observations we might have,

$\begin{array}{c}{p}_N\left(p\text{|}x\right)=\prod _{q=1}^N\rho \left(p\text{|}x_q\right)\\ =\rho \left(p\text{|}x_N\right)p_{N-1}\left(p\text{|}x\right).\end{array}$

We then can exponentially devalue stale information by updating in the log-domain,

log(p_N(p|x))=log(ρ(p|x_N))+λ log((p_N−1(p|x))+c,

where 0<λ<1 and c is a normalization constant so that p_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 _N =p _N−1+v.

The model then seeks to learn v and to update p_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_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 _off =u+e ₀

x _on =u+h+e ₁,

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₁ is the noise vector for the new signal present case, e₀ 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˜N (ū,R_uu) and h˜N (h,R_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_hh+R_ee)N(u−ū,R_uu)

p(x,u|OFF)=N(x−u,R_ee)N(u−ū,R_uu).

After integrating out u one obtains,

p(x|ON)=N(x−h−ū,R_uu+R_hh+R_ee)

p(x|OFF)=N(x−ū,R_uu+R_ee)

From this we can see that the posterior probability of x adhering to the ON model is,

$p\left(\mathrm{ON}x\right)=\frac{N\left(x-\stackrel{\_}{h}-\stackrel{\_}{u},R_{\mathrm{uu}}+R_{\mathrm{hh}}+R_{\mathrm{ee}}\right)p_1}{N\left(x-\stackrel{\_}{h}-\stackrel{\_}{u},R_{\mathrm{uu}}+R_{\mathrm{hh}}+R_{\mathrm{ee}}\right)p_1+N\left(x-\stackrel{\_}{u},R_{\mathrm{uu}}+R_{\mathrm{ee}}\right)p_0},$

where p₁ is the prior probability that we see a new signal arriving and p₀=1−p₁. 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(ū₁,R₁₁),

ū ₁ =ū+R _uu(R_uu+R_hh+R_ee)^-1(x−ū−h)

R ₁₁ =R _uu −R _uu(R_uu+R_hh+R_ee)^-1R_uu;

p(u|x,OFF)=N(ū₀,R₀₀),

ū ₀ =ū+R _uu(R_uu+R_ee)^-1(x−ū)

R ₀₀ =R _uu −R _uu(R_uu+R_ee)^-1R_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_q, qϵ 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_q that correspond to a shape category c_e. The identification problem also often has a time dependency, so that one typically makes an identification over a time period tϵ[r−T,τ]. 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,

$p_k=\sum _{j\in \mathrm{par}\left(k\right)}p_{kj}p_j.$

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,

$x_k=f\left(\sum _{j\in \mathrm{par}\left(k\right)}w_{\mathrm{kj}}x_j\right).$

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_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,

$f\left(x_1,x_2,\dots ,x_I\right)\approx \sum _q{\sigma }_q\prod _if_{\mathrm{qi}}\left(x_i\right).$

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_0,3, the microscopic Maxwell's equations can be written as,

$\begin{array}{cc}\left(\frac{1}{c}{\partial }_t+D\right)F=\frac{1}{c{\varepsilon }_0}\left(c\rho -j\right),& \left(14\right)\end{array}$

where c is the velocity of light, ∂_t is the partial derivative operator with respect to time, D≡e₁∂_x+e₂∂_y+e₃∂_z is the Dirac operator, F≡E+icB is the electromagnetic field vector (or Faraday vector), ρ is the charge density and j is the current density vector. A vector is expressed by v≡Σ_k=1³v^ke_k, and the coordinate vectors e_k have an associative product that has the properties e_qe_m=−e_me_q for q≠m and e_q²=1. The pseudoscaler i≡e₁e₂e₃ commutes with the e_k and has the structural property i²=−1.

A product of vectors in the geometric algebra includes both the inner product and the cross product via the formula,

ab=a·b+ia×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πjft), we see that the electromagnetic field vector satisfies the Helmholtz equation,

DF=−jkF,

where the wavenumber

$k=\frac{1}{c}2\pi f.$

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²=∂_x²+∂_y²+∂_z² is the scalar Laplace operator of the standard wave equation.

The Cauchy integral formula for our Geometric algebra can be written in general by,

$\begin{array}{cc}{\int }_{\partial R}{\mathrm{hdx}}^{*}g={\int }_R\left(\left(\mathrm{hD}\right)g+\mathrm{hDg}\right)\mathrm{dV},& \left(15\right)\end{array}$

where dx*≡dydze₁−dxdze₂+dxdye₃=Σ_kdx*_ke_k, and ⋅* is the Hodge star operator, R is a volume in three space and ∂R is the boundary of the volume. The vector form dx* will be normal to any parameterized surface. In (15) if we substitute h→b and g→Da, we get,

${\int }_{\partial R}{\mathrm{bdx}}^{*}\mathrm{Da}={\int }_R\left(\left(\mathrm{bD}\right)\mathrm{Da}+{\mathrm{bD}}^2a\right)\mathrm{dV}.$

Similarly if we substitute h→bD and g→a we get,

${\int }_{\partial R}{\mathrm{bDdx}}^{*}a={\int }_R\left(\left(D^2b\right)a+\mathrm{DbDa}\right)\mathrm{dV}.$

Subtracting these two equations yields a second order Green's formula,

${\int }_{\partial R}\left({\mathrm{bDdx}}^{*}a-{\mathrm{bdx}}^{*}\mathrm{Da}\right)={\int }_R\left(\left(D^2b\right)a-{\mathrm{bD}}^2a\right)\mathrm{dV}.$

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,

$\begin{array}{cc}{\int }_R\left(\left(h\left(D-\mathrm{jk}\right)\right)g+h\left(D+\mathrm{jk}\right)g\right)\mathrm{dV}={\int }_{\partial R}{\mathrm{hdx}}^{*}g,& \left(16\right)\\ {\int }_{\partial R}\left(b\left(D-\mathrm{jk}\right){\mathrm{dx}}^{*}a-{\mathrm{bdx}}^{*}\left(D-\mathrm{jk}\right)a\right)=\dots {\int }_R\left(\left(\left(D^2+k^2\right)b\right)a-b\left(D^2+k^2\right)a\right)\mathrm{dV}.& \left(17\right)\end{array}$

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 ₊=½(1+ij), and ₋=½(1−ij). The algebra isomorphism is better understood by noting,

_{+ +}=₊

_{− −}=₋

_{+ −}=0

₊+₋=1

i ₋ −i ₊ =j

₋ e ^iwt+₊e^-iwt=e^jwt.

Thus, we can write the Clifford electromagnetic field as,

$\begin{array}{c}F\left(k,r\right)e^{j\omega t}=\mathrm{\_F}\left(k,r\right)e^{i\omega t}+{}_{+}\stackrel{\_}{F}\left(k,r\right)e^{-i\omega t}\\ ={}_{-}F\left(k,r\right)e^{i\omega t}+{}_{+}F\left(-k,r\right)e^{-i\omega t},\end{array}$

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

$k=\frac{2\pi \omega }{c},\mathrm{and}-k\mathrm{or}\frac{-2\pi \omega }{c}.$

One can recover the original representation by replacing,

₊→½(1+ij)

₋→½(1−ij),

resulting in,

$\begin{array}{c}F\left(k,r\right)e^{j\omega t}=\frac{1}{2}\left(F\left(k,r\right)e^{i\omega t}+\stackrel{\_}{F}\left(k,r\right)e^{-i\omega t}\right)+j\frac{1}{2}(F\left(k,r\right)e^{i\omega t}-\\ \stackrel{\_}{F}\left(k,r\right)e^{-i\omega t})\\ =\left(F_r+{\mathrm{jF}}_i\right)e^{j\omega t},\end{array}\hspace{1em}$

where,

F _r=½(F(k,r)+F(k,r))

F _i=½(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=l and let g=F, the electro-magnetic field vector, so that (D+jk)

$F=\frac{1}{c{\varepsilon }_0}\left(c\rho -j\right).$

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,

$\begin{array}{c}\frac{A}{c{\varepsilon }_2}\left(c{\rho }_s-j_s\right)-{\int }_R\mathrm{jkF}={\int }_{\partial R}{\mathrm{dx}}^{*}F\\ ={\int }_{\partial R}\mathrm{nFdS}\\ =\mathrm{An}\left(F_1-F_2\right),\end{array}\hspace{1em}$

where ρ_s is the surface charge density and j_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→0, we also have, ∫_RjkF→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,

$n\left(F_1-F_2\right)=\frac{1}{{\varepsilon }_2}{\rho }_s-c{\mu }_2j_s.$

This shows in particular that,

n×(E₁−E₂)=0

n·(B₁−B₂)=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₋F≡½(F+n{tilde over (F)}n),

where ã is the principle involution defined by {tilde over (e)}_k=−e_k for k=1 . . . 3, and ãb=ã{tilde over (b)}. Thus the continuity relations can be written as,

Q ₋(F₁−F₂)=0.

If you define the complementary projection operator onto the normal component of E by,

Q ₊ F=½(F−n{tilde over (F)}n),

then we have on the surface,

$\begin{array}{c}\left(Q_{+}+Q_{-}\right)F=F\\ Q_{+}^2=Q_{+}\\ Q_{-}^2=Q_{-}\\ Q_{+}\left(F_1-F_2\right)=\frac{1}{{\varepsilon }_2}{\rho }_sn-c{\mu }_2n\times j_s\\ =n\left(\frac{1}{{\varepsilon }_2}{\rho }_s-c{\mu }_2j_s\right)\end{array}\hspace{1em}$

For a perfect conductor the surface current paravector

$K=\frac{1}{{\varepsilon }_2}{\rho }_s-c{\mu }_2j_s,$

can be found from the incident Faraday vector F₀ by,

K=2nF₀,

where is the real part of the Clifford number a with respect to the pseudoscalar i=e₁e₂e₃. 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 _r(x)=−n{tilde over (F)}₀(n{tilde over (x)}n)n,

Also of interest is the reversion involution defined by (ab)^†=b^†a^† and e_k^†=e_k. Note that ĩ=−i, and i^†=−i. Also a=½(a+a^†). For non-deal conductors we will generalize this rule somewhat via an unknown reflection coefficient α, so that,

F _r(x)=−αn{tilde over (F)}₀(n{tilde over (x)}n)n   (18)

The reflection coefficient α 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₀){tilde over ( )}n+x₀, for some x₀ on the reflecting surface.

This model yields the following surface current paravector (vector plus a scalar),

K=2nF₀α^†  (19)

where F₀ 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,

$\left(\frac{1}{2}{\partial }_t+D\right)F_r\left(x\right)e^{j\omega t}=0,$

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 α is a property of the underlying material. Since α can be learned during calibration and estimated from the received data, changes in α 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₀)(D−jk)=δ(x−x₀).

Furthermore consider a Faraday vector g(x) satisfying the Helmholtz version of Maxwell's equation,

$\left(D+\mathrm{jk}\right)g\left(x\right)=\frac{1}{c{\varepsilon }_0}\left(c\rho -j\right).$

Substitution into (16) yields,

${\int }_R\left(\delta \left(x-x_0\right)g+h\left(x-x_0\right)\frac{1}{c{\varepsilon }_0}\left(c\rho \left(x\right)-j\left(x\right)\right)\right)\mathrm{dV}={\int }_{\partial R}h\left(x-x_0\right){\mathrm{dx}}^{*}g\left(x\right).$

Assuming that x₀ lies in the spatial region R (possibly unbounded assuming decaying radiation conditions), we have,

$g\left(x_0\right)={\int }_{\partial R}h\left(x-x_0\right){\mathrm{dx}}^{*}g\left(x\right)-{\int }_Rh\left(x-x_0\right)\frac{1}{c{\varepsilon }_0}\left(c\rho \left(x\right)-j\left(x\right)\right)\mathrm{dV}.$

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,

$\begin{array}{cc}g\left(x_0\right)={\int }_{\partial R}h\left(x-x_0\right){\mathrm{dx}}^{*}g\left(x\right).& \left(20\right)\end{array}$

The reproducing kernel for geometric algebra can be written as,

$h\left(x\right)=\left(D+\mathrm{jk}\right)\frac{e^{-\mathrm{jkr}}}{-4\pi r},r\equiv \sqrt{x^2+y^2+z^2}={\left(\mathrm{xx}\right)}^{\frac{1}{2}}.$

This can be written as,

$\begin{array}{cc}h\left(x\right)=\frac{-e^{-\mathrm{jkr}}}{4\pi r}\left(\mathrm{jk}-\left(\mathrm{jk}+\frac{1}{r}\right)\frac{x}{r}\right)& \left(21\right)\\ x\equiv {\mathrm{xe}}_1+{\mathrm{ye}}_2+{\mathrm{ze}}_3.& \end{array}$

For the Green's theorem in (17), it is helpful to first define the Electromagnetic paravector potential as,

a=cA−ϕ,

where A is the electro-magnetic vector potential and where ϕ is the scalar potential. We also assume the Lorenz gauge which forces

${〈\left(D-\frac{1}{c}\frac{\partial }{\partial t}\right)a〉}_S=0,$

where ⋅_S is the scalar part of a Clifford number. This means,

$\mathrm{cD}·A+\frac{1}{c}\frac{\partial \phi }{\partial t}=0.$

Maxwell's equations now give us,

$\begin{array}{cc}\left(D-\frac{1}{c}\frac{\partial }{\partial t}\right)a=E+\mathrm{icB}=F.& \left(22\right)\end{array}$

For the time harmonic case we get the Helmholtz version of this as,

$\left(D-\mathrm{jk}\right)a=F\left(D^2+k^2\right)a=\frac{1}{c{\varepsilon }_0}\left(c\rho -j\right).$

In Equation (17) we now use a as the electromagnetic paravector potential and b as a Green's function satisfying (D²+k²)b (x−x₀)=δ(x−x₀). This yields,

$\underset{\partial R}{\int }\left(\left(b\left(x-x_0\right)\left(D-\mathrm{jk}\right)\right){\mathrm{dx}}^{*}a\left(x\right)-b\left(x-x_0\right){\mathrm{dx}}^{*}F\left(x\right)\right)+\cdots \hspace{1em}{\int }_R\frac{b\left(x-x_0\right)}{c{\varepsilon }_0}\left(c\rho \left(x\right)-j\left(x\right)\right)\mathrm{dV}=a\left(x_0\right).$

It is customary to choose a Green's function that is constant or even zero on ∂R, however our typical use case evaluates the paravector potential at an arbitrary point x₀ in space, requiring our surface to vary with the choice of receiver location x₀. For this reason it suffices to simply choose the standard Green's function for the sphere,

$b\left(x-x_0\right)=\frac{-e^{-\mathrm{jkr}}}{4\pi r},r=x-x_0$

from which we can obtain,

$\begin{array}{cc}\mathrm{h\_}\left(x-x_0\right)\equiv b\left(x-x_0\right)\left(D-\mathrm{jk}\right)=\frac{-e^{-\mathrm{jkr}}}{4\pi r}\left(-\mathrm{jk}-\left(\mathrm{jk}+\frac{1}{r}\right)\frac{x-x_0}{r}\right).& \left(23\right)\end{array}$

In a region with no sources we can therefore write the reproducing kernel for the paravector potential as,

$\begin{array}{cc}a\left(x_0\right)=\underset{\partial R}{\int }\left(\mathrm{h\_}\left(x-x_0\right){\mathrm{dx}}^{*}a\left(x\right)-b\left(x-x_0\right){\mathrm{dx}}^{*}F\left(x\right)\right).& \left(24\right)\end{array}$

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,

${\int }_R\frac{b\left(x-x_0\right)}{c{\varepsilon }_0}\left(c\rho \left(x\right)-j\left(x\right)\right)\mathrm{dV}=a\left(x_0\right).$

Also note that the scalar potential can, in particular can be written as,

$\phi \left(x_0\right)=\underset{\partial R}{\int }{〈\mathrm{h\_}\left(x-x_0\right){\mathrm{dx}}^{*}a\left(x\right)-b\left(x-x_0\right){\mathrm{dx}}^{*}F\left(x\right)〉}_S.$

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 ∂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 ∂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,

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 α_qF_r(x_q)=−α_qn_q{tilde over (F)}₀(p+r_q)n_q, where F₀(x) is the transmitted Faraday vector seen at position x, n_q is the unit normal vector seen at position x_q=p+r_qϵ∂R.

3. Compute the Faraday vector at antenna position y using an approximation to the reproducing kernel formula,

$g\left(y_m\right)\approx \sum _qh\left(x_q-y_m\right)n_q{\alpha }_aF_r\left(x_q\right)A_q,$

where A_q is the area associated with the sampled surface point x_q.

4. Compute per antenna voltage h_m≈b_m⋅E(y_m), where b_m is the appropriate vector effective length for receiver antenna m, and E(y_m) is the electric field component of g(y_m). Computations can be reduced by bringing the E field projection and the inner product into the sum/integral, ie

$h_m\approx \sum _q{\alpha }_qA_qb_m·{〈h\left(x_q-y_m\right)n_qF_r\left(x_q\right)〉}_{\mathrm{RV}},$

where u_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_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,

ϑ_mq≡A_qb_m·h(x_q−y_m)n_qF_r(x_q)_RV,

in advance after calibration. In this mode we would compute ϑ_mq over a dense grid x_q, indexed by q. Thus with the ϑ_mq pre computed we can write an approximation for the channel vector as,

$\begin{array}{cc}{h}_m=\sum _{q:x_q\in \partial R}{\vartheta }_{\mathrm{mq}}{\alpha }_q.& \left(25\right)\end{array}$

In this model we might consider the reflection coefficients α_q and even the boundary set ∂R as free parameters, not to mention the implicit dependency on the surface normal unit vectors n_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 following algorithm,

Surface Inference Algorithm

1. Hypothesize a boundary region ∂R, containing a center as part of the parameter vector p.

2. Identify all nodes qϵ∂R and compute or update the unit normal vectors n_q if necessary, based on the shape of ∂R.

3. Compute h_m=Σ_q:xqϵ∂Rϑ_mqα_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 α_q in a locally constrained manner (ie small variance), and use nearby q to determine if the bondary ∂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^-jkr, where r is the total distance traveled for the transmitted waveform (about twice the distance to the reflector). The

$\frac{e^{-\mathrm{jkr}}}{r}$

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₀, which we also assume is the transmitter position for simplicity. Assume also that x_m=x₀+d_m is the position of receive antenna m. The total distance traveled to antenna m is,

r _m =y ₀ +∥x ₀ −y+d _m∥.

The total distance to x₀ is,

$\begin{array}{c}{r}_0=2y_0\\ \equiv 2x_0-y.\end{array}$

If we factor out the phase center we have the relative phase term,

e_rel≡e^-jk(rm-r0).

But now we can write,

$\begin{array}{c}{r}_m-r_0=x_0-y+d_m-y_0\\ =\frac{d_m\left(\frac{2d_m}{d_m}.\left(x_0-y\right)+d_m\right)}{x_0-y+d_m-y_0}.\end{array}$

The last expression avoids the loss of precision inherent in subtracting two large numbers r_m−r₀ 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^-jkr0 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

$\left(x,y,z\right)$ $D=e_1\frac{\partial }{\partial x}+e_2\frac{\partial }{\partial y}+e_3\frac{\partial }{\partial z}.$

This coordinate system is good for representing rectangular surfaces and boxes.

Cylindrical Coordinate System

$x=\mathrm{\rho cos}\left(\theta \right)$ $y=\mathrm{\rho sin}\left(\theta \right)$ $z=z,\rho =\sqrt{x^2+y^2}$ $\theta ={\tan }^{-1}\left(\frac{y}{x}\right)$ $z=z,D=D\rho \frac{\partial }{\partial \rho }+D\rho \frac{\partial }{\partial \theta }+D\rho \frac{\partial }{\partial z}$ $D\rho =\cos \left(\theta \right)e_1+\sin \left(\theta \right)e_2$ $D\theta =\frac{-\sin \left(\theta \right)}{\rho }e_1+\frac{\cos \left(\theta \right)}{\rho }e_2.$

This coordinate system is good for representing cylindrical surfaces.

Spherical Coordinate System

$x=r\sin \left(\theta \right)\cos \left(\phi \right)$ $y=r\sin \left(\theta \right)\sin \left(\phi \right)$ $z=r\cos \left(\theta \right)$ $r=\sqrt{x^2+y^2+z^2}$ $\theta ={\cos }^{-1}\left(\frac{z}{r}\right)$ $\phi ={\tan }^{-1}\left(\frac{y}{x}\right)$ $D=\mathrm{Dr}\frac{\partial }{\partial r}+D\theta \frac{\partial }{\partial \theta }+D\phi \frac{\partial }{\partial \phi }$ $\mathrm{Dr}=\sin \left(\theta \right)\cos \left(\phi \right)e_1+\sin \left(\theta \right)\sin \left(\phi \right)e_2+\cos \left(\theta \right)e_3$ $D\theta =\frac{\cos \left(\theta \right)\cos \left(\phi \right)}{r}e_1+\frac{\cos \left(\theta \right)\cos \left(\phi \right)}{r}e_2-\frac{\sin \left(\theta \right)}{r}e_3$ $D\phi =\frac{-\sin \left(\phi \right)}{r\sin \left(\theta \right)}e_1+\frac{\cos \left(\phi \right)}{r\sin \left(\theta \right)}e_2$

This coordinate system is good for representing spheres.

Prolate Spheroidal Coordinate System

$x=a\sqrt{\left({\sigma }^2-1\right)\left(1-r^2\right)}\cos \left(\phi \right)$ $y=a\sqrt{\left({\sigma }^2-1\right)\left(1-r^2\right)}\sin \left(\phi \right)$ $z=a\sigma \tau ,\sigma =\frac{1}{2a}\left(\sqrt{x^2+y^2+{\left(z+a\right)}^2}+\sqrt{x^2+y^2+{\left(z-a\right)}^2}\right)$ $\tau =\frac{1}{2a}\left(\sqrt{x^2+y^2+{\left(z+a\right)}^2}-\sqrt{x^2+y^2+{\left(z-a\right)}^2}\right)$ $\phi ={\tan }^{-1}\left(\frac{y}{x}\right)$ $D=D\sigma \frac{\partial }{\partial \sigma }+D\tau \frac{\partial }{\partial \tau }+D\phi \frac{\partial }{\partial \phi }$ $D\sigma =\frac{x\left(r_{+}+r_{-}\right)}{2{\mathrm{ar}}_{+}r_{-}}e_1+\frac{y\left(r_{+}+r_{-}\right)}{2{\mathrm{ar}}_{+}r_{-}}e_1+\left(\frac{z\left(r_{+}+r_{-}\right)-a\left(r_{+}-r_{-}\right)}{2{\mathrm{ar}}_{+}r_{-}}\right)e_3$ $D\tau =\frac{-x\left(r_{+}-r_{-}\right)}{2{\mathrm{ar}}_{+}r_{-}}e_1+\frac{-y\left(r_{+}-r_{-}\right)}{2{\mathrm{ar}}_{+}r_{-}}e_1+\left(\frac{-z\left(r_{+}-r_{-}\right)+a\left(r_{+}+r_{-}\right)}{2{\mathrm{ar}}_{+}r_{-}}\right)e_3$ $D\phi =\frac{-\sin \left(\phi \right)}{a\sqrt{\left({\sigma }^2-1\right)\left(1-{\tau }^2\right)}}e_1+\frac{\cos \left(\phi \right)}{a\sqrt{\left({\sigma }^2-1\right)\left(1-{\tau }^2\right)}}e_2,\mathrm{where},r_{+}\equiv \sqrt{x^2+y^2+{\left(z+a\right)}^2}$ $r_{-}\equiv \sqrt{x^2+y^2+{\left(z-a\right)}^2},$

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 ξ₁ is constant,

$\frac{x^2}{a^2-{\xi }_1}+\frac{y^2}{b^2-{\xi }_1}+\frac{z^2}{c^2-{\xi }_1}=1$

hyperboloids of one sheet when ξ₂ is constant,

$\frac{x^2}{a^2-{\xi }_2}+\frac{y^2}{b^2-{\xi }_2}+\frac{z^2}{{\xi }_2-c^2}=1$

and hyperboloids of two sheets when ξ₃ is constant

$\frac{x^2}{a^2-{\xi }_3}+\frac{y^2}{b^2-{\xi }_3}+\frac{z^2}{{\xi }_3-c^2}=1,$

where we additionally require,

ξ₁<c²<ξ₂<b²<ξ₃<α².

Coordinates can be derived as,

$x^2=\frac{\left(a^2-{\xi }_1\right)\left(a^2-{\xi }_2\right)\left(a^2-{\xi }_3\right)}{\left(a^2-b^2\right)\left(a^2-c^2\right)}$ $y^2=\frac{\left(b^2-{\xi }_1\right)\left(b^2-{\xi }_2\right)\left({\xi }_3-b^2\right)}{\left(a^2-b^2\right)\left(b^2-c^2\right)}$ $z^2=\frac{\left(c^2-{\xi }_1\right)\left(c^2-{\xi }_2\right)\left({\xi }_3-c^2\right)}{\left(a^2-c^2\right)\left(b^2-c^2\right)}.$

The Dirichlet operators for the confocal ellipsoidal coordinates must be derived implicitly since (ξ₁,ξ₂,ξ₃) can only be found from the roots of the third order polynomial equation in z,

$\frac{x^2}{a^2-z}+\frac{y^2}{b^2-z}+\frac{z^2}{c^2-z}=1.$

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,

$D=D{\psi }_1\frac{\partial }{\partial {\psi }_1}+D{\psi }_2\frac{\partial }{\partial {\psi }_2}+D{\psi }_3\frac{\partial }{\partial {\psi }_3},$

for a coordinate parameterization (ψ₁,ψ₂, ψ₃) . Suppose we know that given a fixed ψ₃=0, the coordinate system forms a closed surface (e.g. for spherical coordinates if ψ₃=r−1, then the closed surface is a sphere).

Define,

$D_{12}=D{\psi }_1\frac{\partial }{\partial {\psi }_1}+D{\psi }_2\frac{\partial }{\partial {\psi }_2}+\mathrm{jk}.$

Furthermore suppose we have a differential operator μ(ψ₃) that depends on ψ₃ and let ϕ(ψ₁,ψ₂) be the value of F restricted to the surface ψ₃−0. We will make the assumption that,

(D+jk)μ(ψ₃)ϕ(ψ₁,ψ₂)=0   (26)

so that F(ψ₁,ψ₂,ψ₃)=μ(ψ₃)ϕ(ψ₁,ψ₂), satisfies the Helmholtz equation. If we require the operator μ(ψ₃) to be valid for any surface function ϕ(ψ₁,ψ₂), then Equation (26) requires the following functional equation,

$\begin{array}{cc}\begin{array}{c}\left(D_{12}+D{\psi }_3\frac{\partial }{\partial {\psi }_3}\right)\mu \left({\psi }_3\right)=0\\ D{\psi }_3{\mu }^{\prime }\left({\psi }_3\right)=-D_{12}\mu \left({\psi }_3\right)\\ {\mu }^{\prime }\left({\psi }_3\right)=-{\left(D{\psi }_3\right)}^{-1}D_{12}\mu \left({\psi }_3\right).\end{array}\hspace{1em}& \left(27\right)\end{array}$

This yields the operator solution to μ(ψ₃),

μ(ψ₃)=exp(−Dψ₃)^-1D₁₂ψ₃).

Thus the solution to the Helmholtz equation becomes,

F(ψ₁,ψ₂,ψ₃)=exp(−(Dψ₃)^-1D₁₂ψ₃)ϕ(ψ₁,ψ₂).

As long as ϕ(ψ₁,ψ₂) is C^∞, we can expand the operator function exp(−(Dψ₃)^-1D₁₂ψ₃) in power series of functions and derivatives with respect to ψ₁ and ψ₂ 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,

$\begin{array}{cc}\mu \left(s\right)=-{\int }_0^s{\left(D{\psi }_3\right)}^{-1}D_{12}\mu \left({\psi }_3\right)d{\psi }_3+1,& \left(28\right)\end{array}$

noting that the initial conditions require μ(0)=1. If D_c≡(Dψ₃)^-1D₁₂, then we can write this as the integro-differential operator,

(1+∫D_c)μ=1,

where ∫ g≡∫₀^s g(ψ₃)dψ₃. One can now argue that,

$\mu =\sum _{q=0}^{\infty }{\left(-\int D_c\right)}^q,$

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

ϕ(ψ₁,ψ₂)=−αn{tilde over (F)}₀(n{tilde over (x)}n)n,

for xϵ∂R defined by ψ₃=0, and where F₀(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_r, we can recover the outgoing surface component of this wave for ψ₃≥0 via,

$P_{+}F_r\equiv \underset{{\psi }_3\to 0_{+}}{\lim }F_r.$

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

μP₊F_r=F_r.

We now involve the surface projection operators, Q₊+Q₋=1.

$\begin{array}{c}\mu \left(Q_{+}+Q_{-}\right)P_{+}F_r=F_r\\ \left(1-\mu Q_{+}P_{+}\right)F_r=Q_{-}P_{+}F_r\\ \left(1-\mu Q_{+}P_{+}\right)F_r=-Q_{-}P_{+}F_0\\ F_r=-{\left(1-\mu Q_{+}P_{+}\right)}^{-1}Q_{-}P_{+}F_0\\ =-\sum _{q=0}^{\infty }{\left(\mu Q_{+}P_{+}\right)}^qQ_{-}P_{+}F_0\end{array}\hspace{1em}$

We used here the electromagnetic continuity condition at the conductor boundary which requires Q₋P₊(F_r+F₀)=0, where F₀ is the incident wave.

This can be combined with Equation (28) to obtain the recursive formulas,

$F_r=\mu \left({\psi }_3\right)\left(Q_{+}P_{+}F_r-Q_{-}P_{+}F_0\right)$ $\mu \left(s\right)=1-{\int }_0^sD_c\mu \left({\psi }_3\right)d{\psi }_3.$

Convergence can be guaranteed for “small” μQ₊P₊, which is possible as ψ₃→0, for sufficiently smooth F₀.

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²+k²)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,

$D^2\psi =\frac{1}{r^2\sin \left(\theta \right)}\left(\sin \left(\theta \right)\frac{\partial }{\partial r}\left(r^2\frac{\partial \psi }{\partial r}\right)+\frac{\partial }{\partial \theta }\left(\sin \left(\theta \right)\frac{\partial \psi }{\partial \theta }\right)+\frac{1}{\sin \left(\theta \right)}\frac{{\partial }^2\psi }{\partial {\phi }^2}\right).$

If we set ψ=R(kr)Θ(θ)ψ(ϕ), into the Helmholtz equation we obtain three separate ordinary differential equations:

$\begin{array}{c}\frac{1}{r^2}\frac{d}{\mathrm{dr}}\left(r^2\frac{\mathrm{dR}}{\mathrm{dr}}\right)+k^2R-\frac{n\left(n+1\right)}{r^2}R=0\\ \frac{1}{\sin \left(\theta \right)}\frac{d}{d\theta }\left(\sin \left(\theta \right)\frac{d\Theta }{d\theta }\right)+\left(n\left(n+1\right)-\frac{m^2}{{\sin }^2\left(\theta \right)}\right)\Theta =0\\ \frac{d^2\Phi }{d\phi }+m^2\Phi =0.\end{array}\hspace{1em}$

From this we can obtain the general solutions,

R(r)=B₁h_n¹(kr)+B₂h_n²(kr),

where h_n¹(kr) is the spherical Hankel function of the first kind, representing an inward traveling wave and h_n²(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,

$h_0^1\left(\mathrm{kr}\right)=-i\frac{\exp \left(\mathrm{ikr}\right)}{\mathrm{kr}}$ $h_0^2\left(\mathrm{kr}\right)=i\frac{\exp \left(-\mathrm{ikr}\right)}{\mathrm{kr}}$ $h_n^1\left(z\right)={z^n\left(-\frac{1}{z}\frac{d}{\mathrm{dz}}\right)}^nh_0^1\left(z\right)$ $h_n^2\left(z\right)={z^n\left(-\frac{1}{z}\frac{d}{\mathrm{dz}}\right)}^nh_0^2\left(z\right).$

The polar angle dependency is solved by,

Θ(θ)=P_n^m(cos(θ)),

where P_n^m(x) is the Legendre polynomial of the first kind.

Φ(ϕ)=exp(mjϕ).

The angular functions are typically collated into a single set of functions known as spherical harmonics Y_n^m(θ,ϕ),

$Y_n^m\left(\theta ,\phi \right)=\sqrt{\frac{\left(2n+1\right)\left(n-m\right)!}{4\pi \left(n+m\right)!}}P_n^m\left(\cos \left(\theta \right)\right)e^{\mathrm{jm}\phi }.$

The spherical harmonics are orthogonal satisfying,

${\int }_0^{2\pi }d\phi {\int }_0^{\pi }Y_{n^{\prime }}^{m^{\prime }}\left(\theta ,\phi \right)Y_n^{m*}\left(\theta ,\phi \right)\sin \left(\theta \right)d\theta ={\delta }_{n^{\prime }n}{\delta }_{m^{\prime }m},$

where δ_m′m is the dirac-delta function. These functions are also complete, allowing the representation of any function on the sphere, f(θ,ϕ) as a series of spherical harmonics. The spherical harmonics are often written as a function of the unit normal vector

$\hat{x}=\frac{x}{x},$

as Y_n^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₀¹(kr) and h₀²(kr) into spherical waves when

$r=x-x_0.\mathrm{Let}x_{>}=\max \left\{x,x_0\right\},x_{>}=\max \left\{x,x_0\right\},\hat{x}=\frac{x}{x},{\hat{x}}_0=\frac{x_0}{x_0}.$

We can then write,

$\begin{array}{cc}\frac{\exp \left(\mathrm{jkr}\right)}{\mathrm{jkr}}=\sum _{l=0}^{\infty }\left(2l+1\right){}_l\left({\mathrm{kx}}_{<}\right)h^1\left({\mathrm{kx}}_{>}\right)P_l\left(\hat{x}·{\hat{x}}_0\right),& \left(29\right)\end{array}$

where j_l(z) is the spherical Bessel function of the first kind and the real part of h¹(z) for real z.

Note that,

r=√{square root over (r₁²+r₀²−2r₁r₀ cos(θ))},

where r₀=∥x₀∥, r₁=∥x∥, and x₀^Tx=r₁r₀ cos(θ).

Suppose that r(r₀,r₁,cos(θ))=√{square root over (r₁²+r₀²−2r₁r₀cos(θ))}, and that we seek to find a decomposition similar to Equation (29) of the form,

$\begin{array}{cc}r\left(r_0,r_1,\cos \left(\theta \right)\right)=\sum _l{\sigma }_l{\alpha }_l\left(r_0\right){\beta }_l\left(r_1\right){\gamma }_l\left(\cos \left(\theta \right)\right).& \left(30\right)\end{array}$

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”, α_l(r₀), β_l(r₁), and γ_l(cos(θ)) . 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,

${}_{j_1j_2\dots j_N}=\sum _{i_1}\sum _{i_2}\dots \sum _{i_N}s_{i_1i_2\dots i_N}U_{j_1i_1}^{\left(1\right)}U_{j_2i_2}^{\left(2\right)}\dots U_{j_Ni_N}^{\left(N\right)},$

where the U_ji⁽ⁿ⁾ are unitary matrices (entry (j,i)), and where s_{i1i2 . . . iN} have all-orthogonal subtensors and ∥s_{i1i2 . . . iN}∥_ij=q² decreases as q increases from 1 to I_j. This means s_{i1i2 . . . iN} is treated as a flat vector while holding one of the indices i_j fixed. The entries s_{i1i2 . . . iN} can be sorted from largest to smallest norm, and the indices flattened to make the decomposition look like,

$\begin{array}{cc}{}_{j_1j_2\dots j_N}=\sum _qs_qU_{j_1i_1\left(q\right)}^{\left(1\right)}U_{j_2i_2\left(q\right)}^{\left(2\right)}\dots U_{j_Ni_N\left(q\right)}^{\left(N\right)}.& \left(31\right)\end{array}$

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₁ to range from 2.2 to 5 meters in increments of a centimeter, and r₀ ranges from 0.2 to 2.0 meters, while cos(θ) 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₁, r₀ and cos(θ) .

Furthermore those functions of r₁, r₀ and cos(θ) 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),

$\begin{array}{cc}{h}_{-}\left(x-x_0\right)=\sum _ld_la_l\left(x\right)b_l\left(y\right)c_l\left(z\right),& \left(32\right)\end{array}$

where x=xe₁+ye₂+ze₃, is the receiver location and d_l are Clifford numbers. The field values are evaluated at the receiver location x₀, 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

$g\left(x_0\right)=\sum _ld_l{\int }_{\partial R}a_l\left(x\right)b_l\left(y\right)c_l\left(z\right){\mathrm{dx}}^{*}F_r\left(x\right).$

The scalar voltage could then be written as,

$V=\sum _lb·{\int }_{\partial R}a_l\left(x\right)b_l\left(y\right)c_l\left(z\right){〈d_l{\mathrm{dx}}^{*}F_r\left(x\right)〉}_{},$

where b is the vector effective length for the receive antenna we are analyzing and F_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 ψ₁,ψ₂,ψ₃, and the closed surface of the scatterer is described by ψ₁=ρ, a constant we can use the HOSVD to decompose the reflected wave

$v\left({\psi }_1,{\psi }_2,{\psi }_3\right)\equiv b·〈{\mathrm{dx}}^{*}F_r\left({\psi }_1,{\psi }_2,{\psi }_3\right)〉,v\left({\psi }_1,{\psi }_2,{\psi }_3\right)\approx \sum _l{\vartheta }_l{\alpha }_l\left({\psi }_l\right){\beta }_l\left({\psi }_2\right){\gamma }_l\left({\psi }_3\right).$

The HOSVD operation is performed over a discretized function of several variables, but we can use splines to interpolate the function values α_l(ψ₁), β_l(ψ₂) and γ_l(ψ₃). We also modify the algorithm to satisfy the Helmholtz equation,

(D²+k²)v(ψ₁,ψ₂,ψ₃)=0

where

$D=D_{{\psi }_1}\frac{\partial }{\partial {\psi }_1}+D{\psi }_2\frac{\partial }{\partial {\psi }_2}+D{\psi }_3\frac{\partial }{\partial {\psi }_3}.$

The operator D² will be a scalar operator that can be discretized and treated as a multi-linear tensor contraction on v(ψ₁,ψ₂,ψ₃). This can be flattened to simply look like a matrix. To obtain the approximation we attempt to solve the Lagrangian problem,

$\begin{array}{cc}\underset{\lambda }{\max }\underset{{\alpha }_l,{\beta }_l,{\gamma }_l,{\vartheta }_l}{\min }{\int }_{{\psi }_1{\psi }_2{\psi }_3}\left({v\left(\rho ,{\psi }_2,{\psi }_3\right)-\sum _l{\vartheta }_l{\alpha }_l\left(\rho \right){\beta }_l\left({\psi }_2\right){\gamma }_l\left({\psi }_3\right)}^2+\dots \lambda {\left(D^2+k^2\right)\sum _l{\vartheta }_l{\alpha }_l\left({\psi }_1\right){\beta }_l\left({\psi }_2\right){\gamma }_l\left({\psi }_3\right)}^2\right)& \left(33\right)\end{array}$

Here the integration ∫_ψ1ψ2ψ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(ψ) is dependent on the center of the scattering body x_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 Σ_lϑ_lα_l(ψ₁)β_l(ψ₂)γ_l(ψ₃) also dependent on x_r. Thus the voltage seen at receiver position x₀ is given by V₀(x_r)=Σ_lϑ_lα_l(ψ₁⁰)β_l(ψ₂⁰)γ_l(ψ₃⁰), where (ψ₁⁰,ψ₂⁰,ψ₃⁰) is the parameterized coordinates of x₀ relative to x_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_r=(x_r,y_r,z_r). This function can be decomposed again using HOSVD techniques to obtain,

$\begin{array}{cc}{V}_0\left(x_r\right)=\sum _l{ϛ}_lq_l\left(x_r\right)u_l\left(y_r\right){\upsilon }_l\left(z_r\right).& \left(34\right)\end{array}$

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,

$\begin{array}{cc}B=\frac{-e^{-\mathrm{jkr}}}{4\pi \mathrm{cr}}\left(\mathrm{kj}+\frac{1}{r}\right)e_r\times b,& \left(35\right)\\ E=\frac{-e^{-\mathrm{jkr}}}{4\pi r}\left(\left(k^2-\frac{\mathrm{jk}}{r}-\frac{1}{r^2}\right)P_{\perp }\left(e_r\right)b+2\left(\frac{\mathrm{kj}}{r}+\frac{1}{r^2}\right)P\left(e_r\right)b\right),& \end{array}$

where e_r≡Dr=x/r and P (e_r)b=e_r(e_r⋅b)=½(b+e_rbe_r) and P⊥(e_r)b=b−e_r(e_r⋅b)=½(b−e_rbe_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₁ 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₂ is it's reflection through the conducting plane P, ie x₂=−nx₁n. Note that the reflection operator is a homomorphism, ie if R_n(x)≡n{tilde over (x)}n, then R_n(ab)=R_n(a)R_n(b), and R_n(x)=x, when xϵP. We can write the reflected field as,

$\begin{array}{c}{F}_r=-n\tilde{F}\left(n\tilde{x}n\right)n\\ =-R_n\left(E\right)+{\mathrm{icR}}_n\left(B\right).\end{array}\hspace{1em}$

Suppose we decompose the E and B field into it's tangential and normal compnueMs,

E=E _n +E _t,

B=B _n +B _t,

then we have,

R _n(E)=−E_n+E_t

R _n(B)=−B_n+B_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,

$\begin{array}{c}{F}_{\mathrm{tot}}=F_r+F\\ =2E_n+2{\mathrm{icB}}_t.\end{array}\hspace{1em}$

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₂=nbn at the reflected position x₂=−nx₁n. The reflected electromagnetic field becomes,

$B=\frac{-e^{-{\mathrm{jkr}}_2}}{4\pi {\mathrm{cr}}_2}\left(\mathrm{kj}+\frac{1}{r_2}\right)e_{r_2}\times b_2,E=\frac{-e^{-{\mathrm{jkr}}_2}}{4\pi r_2}\left(\left(k^2-\frac{\mathrm{jk}}{r_2}-\frac{1}{r_2^2}\right)P_{\perp }\left(e_{r_2}\right)b_2+2\left(\frac{\mathrm{kj}}{r_2}+\frac{1}{r_2^2}\right)P\left(e_{r_2}\right)b_2\right),\mathrm{where},r_2={x-x_2}^2$ $e_{r_2}=\frac{x-x_2}{r_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_n(F).

Conformal Mappings

A conformal map in a Clifford algebra takes the form of a Möbius transform,

ϕ(x)=(ax+b)(cx+d)^-1.

Where a, b, c and d are members of the Clifford group Γ₃, which are products of non-zero vectors in Cl_0,3. These mappings preserve the surface integrals in (16), namely

${\int }_{\partial R}h\left(y\right){\mathrm{dy}}^{*}g\left(y\right)={\int }_{\partial R}h\left(\phi \left(x\right)\right)J^{†}\left(\phi ,x\right){\mathrm{dx}}^{*}J\left(\phi ,x\right)g\left(\phi \left(x\right)\right),$

where

$J\left(\phi ,x\right)\equiv \frac{{\left(\mathrm{cx}+d\right)}^{†}}{{\mathrm{cx}+d}^3}.$

This implies,

${\int }_R\left(\left(h\left(D-\mathrm{jk}\right)\right)g+h\left(D+\mathrm{jk}\right)g\right)\mathrm{dV}={\int }_{{\phi }^{-1}\left(\partial R\right)}h\left(\phi \left(x\right)\right)J^{†}\left(\phi ,x\right){\mathrm{dx}}^{*}J\left(\phi ,x\right)g\left(\phi \left(x\right)\right).$

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,

$g\left(y_0\right)={\int }_{{\phi }^{-1}\left(\partial R\right)}h\left(\phi \left(x\right)-\phi \left(x_0\right)\right)J^{†}\left(\phi ,x\right){\mathrm{dx}}^{*}J\left(\phi ,x\right)g\left(\phi \left(x\right)\right).$

where y₀=ϕ(x₀). This demonstrates that one can use conformal maps to change the surfaces over which a solution can be evaluated.

In particular the Möbius 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 (ϕ,x) g (ϕ(x))=0. However we can use the Vekua operator to map this to a solution that also satisfies the Helmholtz wave equation D²u(x)+k²u(x)=0, via,

$u\left(x\right)=J\left(\phi ,x\right)g\left(\phi \left(x\right)\right)-{\int }_0^1\frac{kx}{2}\frac{t^{\frac{N-2}{2}}}{\sqrt{1-t}}J_1\left(kx\sqrt{1-t}\right)g\left(\phi \left(\mathrm{tx}\right)\right)\mathrm{dt},$

where J₁(⋅) 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_mk^q(p_q) be the voltage seen by the m′th sensor, k′th frequency/channel and the q′th calibration collect, and where p_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,

$\begin{array}{cc}\mu =\sum _{\mathrm{mkq}}{{\upsilon }_{\mathrm{mk}}^q\left(p_q\right)-\sum _l{\alpha }_{\mathrm{mk}}^lb_{lm}\left(p_q\right)}^2.& \left(36\right)\end{array}$

The basis functions b_lm(p_q) are a function of the position/shape parameter vector p_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,

$\sum _qb_{lm}^{*}\left(p_q\right)\left(\sum _n{\alpha }_{\mathrm{mk}}^nb_{nm}\left(p_q\right)-{\upsilon }_{\mathrm{mk}}^q\left(p_q\right)\right)=0$ $\sum _n\sum _qb_{lm}^{*}\left(p_q\right)b_{nm}\left(p_q\right){\alpha }_{\mathrm{mk}}^n=\sum _qb_{lm}^{*}\left(p_q\right){\upsilon }_{\mathrm{mk}}^q\left(p_q\right).$

So for each antenna-frequency pair we have the following L×L linear equation solution,

a_mk=R_bmbm^-1R_bmvmk   (37)

where the l′th element of a_mk is α_mk^l the l, n′th element of the matrix R_bmbm is Σ_qb*_lm (p_q)b_nm (p_q) and the l′th element of R_bmvmk is Σ_qb*_lm (p_q) v_mk^q (p_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 _u(x)=h_n(u)²(kr)Y_n(u)^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_n(u)¹ (kr) for incoming waves. Similar considerations can be made for other coordinate systems.

From this we obtain the electric field from Equation (22) as,

$E_u\left(x\right)=-D{\phi }_u\left(x\right)-{\mathrm{jkcA}}_u\left(x\right)$ ${\phi }_u\left(x\right)=\frac{-c}{\mathrm{jk}}D·A_u\left(x\right)$ $B_u\left(x\right)=D\times A_u\left(x\right).$

We then use one of the several scattering modeling techniques described in Analytic Scattering Models to obtain an induced scattered field Ĕ_u(x_m;p_q) at receiver position x_m and scatterer parameter vector p_q. The scattered field itself may require a series solution for each transmitted E_u(x) and we may see a multi-parameter basis function of the form, Ĕ_pu(x_m;p_q), whose indices (p,u) can either be index flattened to Ĕ_l(x_m; p_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 _lm(p_q)=b_m·Ĕ_l(x_m;p_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:

α_mk^l=β_m^lγ_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 β_m^l terms constant and minimizing over γ_k, followed by holding the γ_k constant and optimizing over β_m^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 }_{\mathrm{mk}}^q\left(p_q\right)\approx \sum _{m^{\prime }}C_{mm^{\prime }}\sum _l{\alpha }_{m^{\prime }k}^lb_{lm^{\prime }}\left(p_q\right).$

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_mk^q(p_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 β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 FFT processing;

the target classification or position parameters are determined from the likelihood function that assumed the scale parameters βq=β, 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 RF IC (ASIC) 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 α, 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 βq over each channel observation q, to absorb bulk phase and gain changes between observations;

the method further includes computing separate channel vectors α_c for each target type c;

the method further includes jointly computing the phase constants βq and the channel vectors α_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 βq=β, 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.

Claims

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 the free parameters of a mathematical model of a channel either as the channel model parameters or in the form of tabulated data;calibrating during the operation of the digital receiver to determine the 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; andsaving target reflector/emitter parameters to be reused for dynamic tracking; and using Bayesian particle filtering or Maximum Likelihood Methods to update mixture models for the 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 the 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 β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 the 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 1, wherein the target classification is determined from the likelihood function after blindly estimating any position dependent parameters during the classification process itself.

11. The method of parameter estimation of claim 9, wherein the target classification is determined from the likelihood function after blindly estimating phase terms such as time delays of arrival.

12. The method of parameter estimation of claim 10, wherein the target classification is determined from the likelihood function after blindly estimating phase delays using FFT processing.

13. The method of parameter estimation of claim 10, wherein the target classification or position parameters are determined from the likelihood function that assumed the scale parameters βq=β, are independent of the 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 the DSP processing hardware is configured to enable parallel likelihood function calculations wherein the same instruction across all processing resources.

15. The method of parameter estimation of claim 1, wherein 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.

16. The method of parameter estimation of claim 1, wherein multiple transmit and receive RF chains are packed into an RF IC (ASIC) 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 1, further interpolating the wave function using basis functions that are reflected through objects of predetermined shape or presumed electrical properties.

19. The method of parameter estimation of claim 17, wherein interpolating uses a likelihood function to compute a fixed set of channel vectors α, modulating the basis functions, that can be used to interpolate the channel at an arbitrary target position.

20. The method of parameter estimation of claim 17, further includes blindly calculating an unknown gain and phase constant βq over each channel observation q, to absorb bulk phase and gain changes between observations.

21. The method of parameter estimation of claim 18, further including computing separate channel vectors αc for each target type c.

22. The method of parameter estimation of claim 18, further including jointly computing the phase constants βq and the channel vectors αc by alternating directions optimization; ie iteratively estimating each parameter while holding the other fixed.

23. The method of parameter estimation of claim 18, 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.

24. The method of parameter estimation of claim 18, further including assuming the scale parameters βq=β, are independent of the collection index, and can therefore be solved using the generalized eigenvalue problem in Equation (8).

25. The method of parameter estimation of claim 18, further calibrating a blind classifier by blindly estimating any position dependent parameters during the classification process itself, without necessarily using known target or array transceiver positions.

26. The method of parameter estimation of claim 18, 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.

27. The method of parameter estimation of claim 18, further calibrating 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.

28. The method of parameter estimation of claim 17, wherein the DSP processing hardware is configured to enable parallel likelihood function calculations, yet can share the same instruction across all processing resources.

29. The method of parameter estimation of claim 17, wherein 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,

30. The method of parameter estimation of claim 17, wherein multiple transmit and receive RF chains are packed into an RF IC, with clocks controlled by a phased locked loop.

31. 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.

32. The method of parameter estimation of claim 17, wherein the likelihood function is itself used to further refine the presumed positions of the transmitter and receiver to improve the calibration channel model.