Cross Reference to Related Inventions

[001] The present invention claims priority to US Provisional Application Serial Number 62/764,814 filed on August 16, 2018, which is incorporated in its entirety by reference.

Background of the Invention

[002] 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.

[003] 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. [004] 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.

[005] 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 herein will apply to both the near and far field and be more accurate due to the removal of information lossy steps.

[006] 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

[007] 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 Figure 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 Figure 2.

[008] The baseline transceiver block diagram is shown in Figure 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.

[009] 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.

[010] Both the target type (shape or dielectric constants etc.) are determined by com puting 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 Figure 4. Each target type must go through an independent calibration procedure, since the target type will determine the basis function coefficients.

[011] Numerous other advantages and features of the invention will become readily ap- parent 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

[012] 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:

[013] FIG. 1 is a radar transceiver;

[014] FIG. 2 is a GPS receiver from multiple satellite transmitters;

[015] FIG. 3 is a transceiver block diagram;

[016] FIG. 4 is a calibration procedure;

[017] FIG. 5 is a delay impulse;

[018] FIG. 6 is a Time delay filtering;

[019] FIG. 7 is a hand wave gesture detection algorithm;

[020] FIG. 8 is a multi-layer neural network;

[021] FIG. 9 is a learning network based on higher order SVDs (HOSVD);

[022] FIG. 10 is a integration box for boundary value calculations;

[023] FIG. 11 shows the HOSVD singular values of a Green’s function;

[024] FIG. 12 is a locus plot of maximal singular functions; and

[025] FIG. 13 shows a dipole image through a reflector.

Description of the Invention

[026] 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.

[027] System Model

[028] 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.

[029] The received signal by the k'th. receiver (or during the k’th observation) is modeled by,

where received waveform vector at sensor k and time is an additive interference plus noise process for sensor is a vector channel process that modulates the transmitted signal 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 and a global parameter vector which is of primary interest to our receiver processing algorithms. The

parameter might consist of geolocation coordinates, target velocity, X-ray mass absorption coefficients, or other signal invariants of interest.

[030] 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, where is the probability of observing the total view column vector

given is the probability density function for the interference

process is the prior distribution of the noisy, free parameters of the channel process

[031] Linear Channel Model

[032] 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,

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

where n is the time sample index, channel matrix as a function of the nuisance parameter vector and parameter of interest vector

for all Q emitters /reflectors. Also the emitter signal waveforms are

packed into the Q x 1 signal vector

[033] Colored Gaussian Interference

[034] Some processing simplifications can be made by subsuming all unmodeled inter fering waveforms into the single interference vector and further presuming that this interference is complex Gaussian, yet not white. In particular it can be presumed the inter- ference takes the form of a zero mean multivariate, complex circularly symmetric Gaussian vector,

The zero mean Gaussian normal distribution is referred to by, and one can also

define the complex normal distribution function as,

This allows the likelihood function to be written as,

where is the matrix trace operator and is the outer product operator.

[035] We now perform some manipulations to expose the sufficient statistics for the prob- ability function as follows:

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

and where A matrix version of completing the square provides a way to identify sufficient statistics for this problem. We write:

Note that is a least squares approximation of the channel matrix, and that is an estimate of the interference covariance matrix. This

allows one to write,

This makes a sufficient statistic given

[036] We can write this in the whitened form, where is the interference whitened and signal power normalized chan

nel matrix for view k, and where is an estimate of the interference

whitened and signal power normalized channel matrix. The matrix is the upper triangu- lar Cholesky factor of being the Hermitian transpose of the inverse of that

matrix. Similarly is the Hermitian transpose of the Cholesky factor of

[037] 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.

[038] This result is dependent on being a sufficiently wide band signal to make

a full rank matrix. We must also cover the very important case of 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 is actually a rank-1 vector modeled as the sum of it’s multipath components,

Here we assume that the position parameters p is shared by all reflected multipath compo- nents.

[039] Narrow Band Channels

[040] The important case wherein 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.

[041] If one specializes to the case wherein k represents a frequency and q is associated with a path delay we can approximate the signal of interest as,

where is a narrow band signal centered at

where Similarly

we have,

[042] This allows Equation (3) to be written as,

where and

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

where

[043] Aperture Noise Model

[044] For proper estimation of channel parameters it is critical to have an accurate noise model of the aperture vector 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.

[045] For simplicity we refer to the parameterized aperture vector as . 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,

From this we can consider several different aperture noise models.

[046] Additive Gaussian Noise [047] This model simply adds Gaussian noise to the true aperture

where e is colored Gaussian complex noise. In general e may be dependent on the unknown parameters It can be useful to consider a piecewise approximation of the e noise pa- rameters 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.

[048] 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.

[049] Multiplicative Aperture Noise

[050] 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, where represents elementwise multiplication and e is colored Gaussian complex noise.

[051] One approximation subsumes this noise model into the simpler noise model of the Additive Gaussian Noise, with e replaced with 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.

[052] 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 vector, where is a possibly whitened component of the channel seen at frequency k across all receivers. The receiver model can be written as,

Here is the whitened received data, in the form of an estimated channel response, playing the role of in Equation (7).

[053] Assuming a complex Gaussian noise model we have,

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

[054] We can now integrate away the nuisance parameter representing the gain to obtain, where

The posterior expectation for and it’s maximum likelihood estimate is given by, It can be shown that to evaluate the argument in the exponential of Equation (9), it suffices to substitute in Equation (8). Also the determinant is easily computed as,

[055] In the case that the multiplicative noise is applied over both antenna number and frequency number, we have the receiver channel vector modeled as, where m is the antenna index and k is the frequency index. This model leads to the conditional distribution of,

where

where is an M x M diagonal matrix, whose m/th diagonal element is and is a diagonal matrix whose m’th diagonal element is the variance of the gain

multiplier

[056] Multi-Target and Multi-Scatter Models

[057] 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),

where 1 is the all ones vector, and where we drop the 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.

[58] From this model we obtain so that we can write the likelihood function as,

[059] 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 /scatter er . 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.

[060] 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.

[061] Multi- Aperture Target and Phase Normalization

[062] 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,

where the target channel is seen as a mixture of multiple related component channels, (e.g. a human body with torso and limbs). The vector is presumed to be the station ary background clutter. For simplicity one can presume that the coefficients are Gaussian distributed, as is the background clutter with,

The probability for the signal on and signal off are,

from this one can compute the posterior probability of the target present (ON) and target not present (OFF), where p ₁ is the prior probability of the target being present and

[063] An additional modification to our model is to allow the removal of bulk phase terms for each channel component. Thus we have, 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.

[064] One can use successive interference cancellation (SIC) or parallel interference cancellation techniques to learn c _q , simply by rewriting,

[065] 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,

We either pick h _p as part of the calibrated array manifold or simply as a raw data estimate.

[066] Position /Parameter Uncertianty

[067] Each channel vector has a dependency on it’s position and other extraneous param eters, 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 may contain the uncertain scatterer / emitter positions. We may assume or estimate priors for these parameters, If these prior distributions are Gaussian, then we have, Defining the real multivariate normal distribution by we can write the likelihood

function as,

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

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.

[069] Range Filtering

[070] 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.

[071] 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, where is the signal delay seen at antenna m, and where is the signal frequency for index k. We perform an inverse discrete Fourier transform, and obtain,

If its assumed that the frequency hopping is linear and that is dominated by a frequency independent mode, we then have,

[072] One can plot what the simple delay channel looks like in the time domain. Here assuming we plot the impulsive response in Figure 5.

[073] As the number of frequencies 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.

[074] The required processing is illustrated in Figure 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.

[075] The time delay filtering concept in Figure 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.

[076] 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. [077] Parameter Estimation and Tracking

[078] 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 esti- mation 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 com putation, either on a multi-core processor or an embedded graphics processing unit (GPU).

[079] Gaussian Mixture Particle Filters

[080] 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.

[081] 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.

[082] Bayes theorem for computing the posterior distribution can be written as,

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

where,

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.

[083] Define the real multivariate Gaussian probability density function as,

One can propose to model the prior and posterior probabilities as a Gaussian mixture. The prior is assumed to adhere to, where One can now propose the following algorithm to update our Gaus

sian mixture so that it approximates the posterior probability given in Equation (12) with

[084] Gaussian Mixture Filter Update

[085] 1. Sample N random points from the Gaussian mixture prior

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

[087] 3. Normalize the likelihoods to sum to one,

[088] 4. Resample a subset of Q < N points from the original N, where the probability of choosing p _n equals Call this set of renumbered points for subindex q and the associated original index,

[089] 5. Use the EM algorithm, and the points to learn a new set of

to fit the set of sampled

[090] We now delve into some of the details of the algorithm.

[091] Sampling the Mixture and Resampling

[092] For a given mixture probability of the form,

we can generate N pseudo-random samples from any such distribution for which we can sample directly from For each n we randomly select what category j, p _n will belong to,

based on the category probabilities Selecting the category can be done in approximately

o(log (J)) comparisons by using a binary search into the cumulative sum array of namely

When a uniform random variate u is chosen, we need merely find j such that,

using binary search, and initializing

[093] Once the category is chosen for a given n, we generate a random sample using the pdf f For example when we can sample from any given multivariate Gaussian distribution by setting where e„ 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 matrix is the upper Cholesky factor of the positive definite covariance matrix R _j . These computations are highly amenable to parallel processing as well

[094] After this we choose to sub-sample our original large sample of N points in propor tion to the likihood

Again we can compute a cumulative sum of the likelihoods followed by

a binary search to find a n where a given, newly generated random uniform variable satisfies, This is our

[095] Let us compute the probability distribution for We have,

We have 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 asymptotically approaches our desired posterior distribution. In what

follows we label the resampled points as

[096] Expectation Maximization Algorithm

[097] The Expectation Maximization (EM) algorithm can be used to learn the free pa rameters 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) .

[098] The Gaussian mixture model for the posterior distribution will take the form,

The initial values for are chosen to be the same as the prior. The EM algo-

rithm consists first performing the expectation of the logarithm of the conditional distribution which is a lower bound to the original function, followed by a maximization of

this lower bound over the unknown parameters to find the next iterate. This

can be summarized as,

[099] Expectation Step

[100] Compute the posterior likelihood of a sampled point p _q being in category j,

The parameter can be thought of as the expected number of points associated to category j. Note

[101] Maximization Step

[102] The updated mixture parameters are given by,

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,

[103] 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 drawn from the prior distribution, modify (13) to read, This causes a given p _n to be directly weighted by it’s likelihood function.

[104] Geolocation and Tracking

[105] 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) .

[106] Dynamic Models

[107] 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 envi- ronment. This can be accommodated by our Bayesian approach, simply by bounding the variance of the parameter of interest from below. So for example if is the posterior distribution of position parameter p given the observed data x, after many observations we might have,

We then can exponentially devalue stale information by updating in the log-domain, where and c is a normalization constant so that 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.

[108] Another approach for system dynamics is to assume that p adheres to a first order discrete constant velocity model, ie

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.

[109] 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.

[110] Detecting Environment Changes

[111] 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,

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, is the signal off noise, and h is the channel vector for a new or moving emitter.

[112] Because we continuously monitor the array we can learn an approximate prior distribution for u and h wherein 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,

After integrating out u one obtains,

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

where p ₁ is the prior probability that we see a new signal arriving and 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.

[113] We also need to compute the posterior distribution for u so that we can update it’s statistics. We have,

The preferred modality is to update only when the OFF state is detected. Note also the sim ilarity of these updates to the Kalman filter. The latter has a close relationship to Bayesian inference, when the distributions in question are normally distributed.

[114] Shape and Gesture Identification

[115] Suppose that a target has a collection of emitters / scatterers associated with it 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 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.

[116] 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. [117] 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 Figure 7.

[118] 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.

[119] 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 Figure 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.

[120] 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.

[121] Using HOSVD in Machine Learning

[122] 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 Figure 8.

[123] 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 prob ability, ie for node k,

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.

[124] 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,

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.

[125] We introduce here a new approach that exploits the structure compression properties of higher order singular value decomposition. In Figure 9 we show that the HOSVD network is essentially a one layer network that learns a set of eigen-functions (the and some multipliers to approximate the nonlinear function of interest. The functions are multiplied and then summed to obtain the final function value via,

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 Figure (11) shows.

[126] Thus there may be significant computational savings in using this flattened struc ture, once the individual eigenfunctions are approximated by splines. This structure provides an innovative way of implementing many of the machine learning networks and for approx- imating the many functions used in partial differential equation boundary value problems, such as the Helmholtz wave equation for electromagnetic propagation.

[127] Channel Modeling

[128] 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 electromagnetic scattering. We propose a series of approximations to the scattering problem that exploit known geometries and some novel expansions of the electromagnetic reproducing kernel.

[129] 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.

[130] Electro-Magnetic Scattering

[131] 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.

[132] 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.

1133] In the geometric Clifford algebra the microscopic Maxwell’s equations can be written as,

where c is the velocity of light, is the partial derivative operator with respect to time,

is the Dirac operator, is the electromagnetic field vector

(or Faraday vector), p is the charge density and j is the current density vector. A vector is expressed by and the coordinate vectors have an associative product

that has the properties and The pseudoscaler

commutes with the e _k and has the structural property

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

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 we see that the electromagnetic field vector satisfies the Helmholtz equation, where the wavenumber Note that the j used here is a square root of -1 for tracking

phase changes for a complex phasor at frequency /, which is different from the pseudoscalar i. An important property of the Dirac operator within a Clifford algebra is the fact that, is the scalar Laplace operator of the standard wave equation.

[135] The Cauchy integral formula for our Geometric algebra can be written in general by, where is the Hodge star operator, R is

a volume in three space and is the boundary of the volume. The vector form dx* will be normal to any parameterized surface. In (15) if we substitute we get,

Similarly if we substitute we get

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

[136] 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,

[137] Bicomplex Algebra

[138] 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.

[139] The two algebras coincide with coefficients of either

. The algebra isomorphism is better understood by noting,

Thus, we can write the Clifford electromagnetic field as,

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 or One can recover the original representation by replacing,

resulting in,

where,

[140] Boundary Conditions

[141] 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 = land let g— F, the electro-magnetic field vector, so that

[142] Let the volume of integration be an infinitely thin box that spans the boundary of two regions as shown in Figure 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,

where is the surface charge density and 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 we also have, 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.

[143] In one fell swoop, this proves the electromagnetic boundary conditions,

This shows in particular that,

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, where a is the principle involution defined by Thus the continuity relations can be written as,

If you define the complementary projection operator onto the normal component of E by, then we have on the surface,

[144] For a perfect conductor the surface current paravector can be found from the incident Faraday vector F ₀ by, where is the real part of the Clifford number a with respect to the pseudoscalar i = Here the fields are evaluated at the surface of the conductor. We can also write

down the reflected field from a perfect conductor as,

Also of interest is the reversion involution defined by Note that

For non-deal conductors we will generalize this

rule somewhat via an unknown reflection coefficient a, so that,

The reflection coefficient a 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 for some x ₀ on the reflecting surface.

[145] This model yields the following surface current paravector (vector plus a scalar), where 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.

[146] It can be shown that the reflected electromagnetic field vector in Equation (18) sat isfies Diracs equation, 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 a is a property of the underlying material. Since a can be learned during calibration and estimated from the received data, changes in a 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.

[147] The Reproducing Kernel

[148] Consider the Helmholtz form of the Cauchy equation in (16). Suppose we can find a function /i(x) in the geometric algebra that has the property that applying the right Helmholtz operator yields the Dirac delta function,

Furthermore consider a Faraday vector g(x) satisfying the Helmholtz version of Maxwell’s equation, Substitution into (16) yields,

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

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,

The reproducing kernel for geometric algebra can be written as,

This can be written as,

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

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

of a Clifford number. This means,

Maxwell’s equations now give us,

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

In Equation (17) we now use a as the electromagnetic paravector potential and b as a Green’s function satisfying This yields,

[150] It is customary to choose a Green’s function that is constant or even zero on 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 XQ. For this reason it suffices to simply choose the standard Green’s function for the sphere, from which we can obtain,

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

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,

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

[151] Reproducing Kernel Channel Model and Surface Inference

[152] 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 dR. 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 dR, 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).

[153] We summarize this algorithm as,

[154] Reproducing Kernel Algorithm [155] 1. Hypothesize a position p for the center of the emitter / reflector region R.

[156] 2. Compute a set of Faraday vectors for the reflected field

where F ₀ (x) is the transmitted Faraday vector seen at position is the unit normal vector seen at position

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

where A _q is the area associated with the sampled surface point x _g .

[158] 4. Compute per antenna voltage where b _m is the appropriate vector effective length for receiver antenna is the electric field component of Computations can be reduced by bringing the E field projection and the inner

product into the sum/integral, ie

where is the real vector part of u. Of course more sophisticated quadrature formulas can be used to obtain the integral approximation.

[159] 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.

[160] Another algorithm suggests itself from the form of the Reproducing Kernel Algo rithm. It is theoretically possible to compute the contribution of all the nodes to the m’th antenna via,

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

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

[161] This suggests the following algorithm,

[162] Surface Inference Algorithm

[163] 1. Hypothesize a boundary region containing a center as part of the parameter vector p.

[164] 2. Identify all nodes and compute or update the unit normal vectors n _q if necessary, based on the shape of

[165] 3. Compute and use it to compute the likelihood for use in

(9) and the algorithm in the Gaussian Mixture Filter.

[166] 4. Use the Gaussian Mixture Filter algorithm to update the mixing gains a _q in a locally constrained manner (ie small variance), and use nearby q to determine if the bondary dR should be updated if it achieves a better fit.

[167] 5. Test convergence if no improvement can be made in the likelihood, if not return to step 1.

[168] Factoring Common Phase Center

[169] For higher frequencies (larger k), spherical harmonics produce phase terms of the form where r is the total distance traveled for the transmitted waveform (about twice the distance to the reflector). The 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.

[170] Suppose the reflector is at position y and the center of the array is at which we also assume is the transmitter position for simplicity. Assume also that is the position of receive antenna m. The total distance traveled to antenna m is,

The total distance to is,

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

But now we can write,

[171] The last expression avoids the loss of precision inherent in subtracting two large numbers 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 in the channel models via relative gain normalization.

[172] Analytic Scattering Models

[173] 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.

[174] We first look at some useful coordinate systems, whose level curves represent objects of interest. We also derive their Dirac operators.

[175] Rectangular Coordinate System

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

[176] Cylindrical Coordinate System

This coordinate system is good for representing cylindrical surfaces.

[177] Spherical Coordinate System

This coordinate system is good for representing spheres.

[178] Prolate Spheroidal Coordinate System

where,

Prolate spheroidal coordinate systems are good for modeling an ellipse that has been rotated about the z-axis.

[179] Confocal Ellipsoidal Coordinates

[180] This confocal ellipsoidal coordinate system has level surfaces corresponding to el lipsoids when is constant,

hyperboloids of one sheet when is constant, and hyperboloids of two sheets when Cs is constant where we additionally require,

Coordinates can be derived as,

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,

[181] Cauchy- Kovalevskaya Extension Theorem

[182] 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, for a coordinate parameterization Suppose we know that given a fixed the coordinate system forms a closed surface (e.g. for spherical coordinates if then the closed surface is a sphere).

[183] Define,

Furthermore suppose we have a differential operator that depends on and let be the value of F restricted to the surface We will make the assumption

that,

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,

This yields the operator solution to Thus the solution to the Helmholtz equation becomes,

As long as we can expand the operator function exp in

power series of functions and derivatives with respect to to obtain a solution that converges locally. This generalizes the Cauchy-Kovalevskaya extension theorem.

[ 184] Another form of this has slightly better convergence properties. Starting from Equation (27) we write, noting that the initial conditions require then we can write this as the integro-difierentiai operator, where One can now argue that,

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 for defined by and where F ₀ (x) is the incident (transmitted) wave.

[185] It is possible to extend the Kovalevskaya extension to actually recover the surface current paravector if 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,

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.

[186] Our extended Kovalevskaya extension recovers the original function from the bound- ary, namely

We now involve the surface projection operators,

We used here the electromagnetic continuity condition at the conductor boundary which requires where F ₀ is the incident wave.

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

Convergence can be guaranteed for“small” which is possible as for suffl- ciently smooth F ₀ .

[188] The main obstacle for use appears to be the proliferation of symbolic terms for larger values of q.

[189] Basis Functions and HOSVD

[190] For a set of basis functions to approximate the scattered field, it is convenient to look the separable solutions of the Helmholtz equation, 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.

[191] The Laplacian in spherical coordinate systems can be written as,

If we set into the Helmholtz equation we obtain three separate

ordinary differential equations:

From this we can obtain the general solutions, where is the spherical Hankel function of the first kind, representing an inward traveling wave and 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,

[192] The polar angle dependency is solved by, where is the Legendre polynomial of the first kind.

[193] The angular functions are typically collated into a single set of functions known as spherical harmonics

The spherical harmonics are orthogonal satisfying,

where is the dirac-delta function. These functions are also complete, allowing the

representation of any function on the sphere, as a series of spherical harmonics. The spherical harmonics are often written as a function of the unit normal vector , as

, which is a polynomial function of the components of , called harmonic polynomials.

[194] It is possible to decompose the spherical green’s functions ) into spherical waves when

We can then write,

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

[195] Note that,

[196] Suppose that and that we seek to find a decomposition similar to Equation (29) of the form,

Furthermore we seek such a decomposition that gives us the best fit to with the min- imum 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”,

This problem can be solved using a device called multilinear sin

gular value decompositions, or higher order svd (HOSVD). This decomposition reduces a multilinear form A (tensor) into a sum of mixed eigenvectors of the form, where the are unitary matrices (entry , and where have all-orthogonal subtensors and decreases as q increases from 1 to . This means is treated as a flat vector while holding one of the indices i _j fixed. The entries can be sorted from largest to smallest norm, and the indices flattened to make the decomposition look like,

The indices for these tensors all represent the values of continuous parameters, that have been densely sampled.

[197] 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 Figure 11. We allow rq to range from 2.2 to 5 meters in increments of a centimeter, and ranges from 0.2 to 2.0 meters, while cos (q) 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 Figure 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 (0) .

[198] Furthermore those functions of r ₁ , r ₀ and cos (q) are smooth as can be seen by the locus plots in Figure 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. [199] 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),

where is the receiver location and 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.

[200] For the reproducing kernel case we have from Equation (20) and Equation (32)for rectangular coordinates

The scalar voltage could then be written as,

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.

[201] 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

The HOSVD operation is performed over a discretized function of several variables, but we can use splines to interpolate the function values We also modify the algorithm to satisfy the Helmholtz equation, where The operator D ² will be a scalar operator that can

be discretized and treated as a multi-linear tensor contraction on This can be flattened to simply look like a matrix. To obtain the approximation we attempt to solve the Lagrangian problem,

Here the integration 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.

[202] Now implicitly is dependent on the center of the scattering body Xr, which is

one of the parameters we need to estimate. After solving the Lagrange constrained minimiza- tion problem in Equation (33) we obtain an approximating series solution

also dependent on x _r . Thus the voltage seen at receiver position x ₀ is given by

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 This function can be decomposed again using HOSVD techniques to obtain,

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.

[203] Tmage Techniques

[204] 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, where,

where

[205] 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, whose superposition with the incident field satisfies

the boundary conditions.

[206] This situation is illustrated in Figure 13.

[207] If the coordinate origin is on the electromagnetic screen, and Xi 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,

Note that the reflection operator is a homomorphism, ie if

when We can write the reflected field as,

Suppose we decompose the E and B field into it’s tangential and normal components, then we have,

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,

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.

[208] 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 at the reflected position The reflected electromagnetic field becomes,

where,

[209] 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

[210] Conformal Mappings

[211] A conformal map in a Clifford algebra takes the form of a Mobius transform,

Where a, b, c and d are members of the Clifford group which are products of non-zero vectors in These mappings preserve the surface integrals in (16), namely where This implies,

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,

where This demonstrates that one can use conformal maps to change the

surfaces over which a solution can be evaluated.

[212] In particular the Mobius transform is known to map spheres to spheres, where a “sphere” of infinite radius is a half plane. Thus we can formulate a solution to the boundary value problem in the half plane, using say the imaging solution of the prior section and then use a conformal map to map the problem to a sphere or vice-versa.

[213] 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

However we can use the Vekua operator to map this to a solution

that also satisfies the Helmholtz wave equation via,

where ) is the first order Bessel function of the first kind. The Vekua transform is actu ally a general way to associate solutions of the ordinary Laplace equation with solutions of the Helmholtz equation.

[214] Calibration

[215] 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.

[216] Let be the voltage seen by sensor, k’th frequency/ channel and

the q’th calibration collect, and where p _q is the“position” and shape parameter vector, for the g’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,

The basis functions are a function of the position/shape parameter vector p _q an

tenna 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. [217] The necessary conditions for the least squares calibration are given by,

So for each antenna-frequency pair we have the following L x L linear equation solution:, where the l’th element of element of the matrix

and the l’th element of

[218] 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.

[219] Obtaining The Basis Functions

[220] 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.

[221] 1. Generate a model for the transmitted electromagnetic field.

[222] 2. Predict the field being scattered (for radar applications) off of the desired target type. [223] 3. Use the prediction to generate basis functions for the received voltage at each antenna receiver and frequency,

[224] 4. Collect data in the near field for multiple scatterer /receiver positions,

[225] 5. Use least squares to estimate the complex combining weights for the basis functions.

[226] 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, 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 for incoming waves. Similar considerations can be made for other coordinate

systems.

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

We then use one of the several scattering modeling techniques described in Analytic Scattering Models to obtain an induced scattered field at receiver position and scatterer parameter vector p _q . The scattered field itself may require a series solution for each transmitted and we may see a multi-parameter basis function of the form,

, whose indices (p, u) can either be index flattened to or further sub-selected on the basis of the largest singular values using HOSVD techniques described in Basis Functions and HOSVD.

[228] Finally the actual voltage and the scalar basis function is obtained via the vector effective length for antenna m,

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.

[229] Calibration Model Modifications

[230] 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.

[231] 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:

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 terms constant and minimizing over followed by holding the constant and optimizing over This is repeated until convergence.

[232] 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 val- ues, 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.

[233] 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,

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

[234] 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.

[235] 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:

[236] the transmitted signal is a frequency stepped radar;

[237] a statistical likelihood function is used to determine target type or target position;

[238] a channel template for each target type is used to further classify the target type or estimate the target position;

[239] static direct path clutter is removed via a cancellation algorithm;

[240] the parameter estimation blindly calculates an unknown gain and phase constant bq over each channel observation q, to absorb bulk phase and gain changes between observations;

[241] 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). [242] the target classification is determined from the likelihood function using the posterior probability of detection via Bayes theorem;

[243] the target classification is determined from the likelihood function after blindly estimating any position dependent parameters during the classification process itself;

[244] the target classification is determined from the likelihood function after blindly estimating phase terms such as time delays of arrival;

[245] the target classification is determined from the likelihood function after blindly estimating phase delays using FFT processing;

[246] the target classification or position parameters are determined from the likelihood function that assumed the scale parameters bq = b, are independent of the collection index and can therefore be solved using the generalized eigenvalue problem in Equation (8);

[247] the DSP processing hardware is configured to enable parallel likelihood function calculations wherein the same instruction across all processing resources;

[248] 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;

[249] multiple transmit and receive RF chains are packed into an RF IC (ASIC) with clocks controlled by a phased locked loop;

[250] 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;

[251] the wave function can be further interpolated using basis functions that are reflected through objects of predetermined shape or presumed electrical properties; [252] interpolating uses a likelihood function to compute a fixed set of channel vectors a, modulating the basis functions, that can be used to interpolate the channel at an arbitrary target position;

[253] the method further includes blindly calculating an unknown gain and phase constant bq over each channel observation q, to absorb bulk phase and gain changes between observations;

[254] the method further includes computing separate channel vectors a _c for each target type c;

[255] the method further includes jointly computing the phase constants bq and the channel vectors a _c by alternating directions optimization; ie iteratively estimating each parameter while holding the other fixed;

[256] 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;

[257] the method further includes assuming the scale parameters bq = b, are independent of the collection index, and can therefore be solved using the generalized eigenvalue problem in

Equation (8);

[258] 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;

[259] 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; [260] 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;

[261] the DSP processing hardware is configured to enable parallel likelihood function calculations, yet can share the same instruction across all processing resources;

[262] 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;

[263] multiple transmit and receive RF chains are packed into an RF IC, with clocks controlled by a phased locked loop;

[264] 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

[265] the likelihood function is itself used to further refine the presumed positions of the transmitter and receiver to improve the calibration channel model.

[266] 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.