Non-linear technique for mitigating correlation timing errors due to multipath signals

[0002] This application claims priority from provisional application No. 60/263,360, filed on Jan. 22, 2001.

[0003] This invention relates to minimization of positioning errors that are experienced by radio locating systems in the presence of multipath signals that arrive from nearby diffractive features and reflective surfaces.

[0004] Radiolocation schemes estimate location of objects to be tracked by measuring time parameters that include Time of Arrival or Time Delay of Arrival at a number of receivers at known locations. Other techniques include measuring Angle of Arrival or a composite of time and angle measures.

[0005] The accuracy of radiolocating system depends on an accurate time measurement of the time of the incident signal sent by the transmitting device to be located by the network of receiving antennas. In real operating environments it is improbable that only the antenna network receives only an incident (line of sight) signal. The transmitted signal also arrives at the antennas from indirect paths. Reflective surfaces, such as buildings, the ground, surrounding vehicles and other surfaces are sources of undesirable signal paths. Also, the transmitted signal may be diffracted by sharp edges of objects. These phenomena result in the transmitted signals traveling to the antenna networks via non-direct paths. The sum of all signal paths from the transmitter to the receivers results in a composite of multipath signals.

[0006] The combined signal at the receivers created by the incident (line of sight) and the diffracted and reflected (non-line of sight) signals creates an uncertainty as to the actual arrival time of the incident signal.

[0007] The present invention overcomes the above-described problems in the prior art by providing a method of determining the time of arrival of an incident signal transmitted utilizing a code word that is interspersed with noise and diffracted and reflected signals of the incident signal. The incident, diffracted and reflected signals and noise are collectively defined as a combined signal. The method comprises receiving the combined signal; obtaining the code word to apply to the combined signal; comparing the combined signal with the code word to filter the incident signal and diffracted and reflected signals arising from the incident signal from the combined signal; utilizing the incident signal and related diffracted and reflected signals to determine a first parameter indicative of the incident signal's time of arrival; and filtering the first parameter utilizing a multipath mitigation algorithm to mitigate the influence of the diffracted and reflected signals to obtain a second parameter indicative of the incident signal's time of arrival.

[0008] The invention will be more readily understood from a reading of the following specifications and by reference to the accompanying drawings forming a part thereof, wherein an example of the invention is shown and wherein:

[0009] FIG. 1 illustrates a typical Transmitter implementation for a Radio Locating System.

[0010] FIG. 2 illustrates a Receiver for a Radiolocation System for mitigating Correlator timing errors in the presence of Multipath Signals.

[0011] FIG. 3 illustrates the Cubic Spline method of estimating a receiver correlation pulse response.

[0012] FIG. 4 illustrates a Normal Autocorrelation Receiver output without received reflections.

[0013] FIG. 5 illustrates a Measured Received Signal and its Cubic Spline Approximation.

[0014] FIG. 6 illustrates a Correlation Noise Peak Model.

[0015] FIG. 7 illustrates the Inphase (I) components.

[0016] FIG. 8 illustrates the equivalent Q components α_Qk·corr(m,m)(t−t_k).

[0017] FIG. 9 compares the original correlator I output (solid) with the mitigation algorithm model (dashed).

[0018] FIG. 10 compares the correlator Q output (solid) with the mitigation algorithm Q model (dashed).

[0019] FIG. 11 compares the correlator I output (solid) with the mitigation algorithm model (dashed).

[0020] FIG. 12 compares the correlator Q output (solid) with the model Q data (dashed).

[0021] FIG. 13 presents both the Inphase (I) components (dotted) and the Quadrature Phase (Q) components (dashed).

[0022] The construction designed to carry out the invention will hereinafter be described, together with other features thereof.

[0023] Turning now to the figures in which like numerals represent like elements throughout the several views, exemplary embodiments of the present invention are described. Radio locating systems use a variety of time measurement techniques that include Time of Arrival, Time Delay of Arrival and Angle of Arrival. This invention will improve the accuracy independently of the time measurement technique.

[0024] FIG. 1 illustrates a typical Transmitter that attaches to a device to be located using radio waves. A source of Information is modulated onto a Carrier that frequency translates the Information to a higher frequency band. Prior to transmission a Pseudo Random Serial Data Stream then modulates the Information at the Carrier Frequency. This modulation process is commonly referred to as Direct Sequence Spread Spectrum (DSSS) or generically, spreading the transmission. A key point is that the common practice is to use a spreading code word that has a high autocorrelation and low cross-correlation characteristic.

[0025] FIG. 1 is a typical implementation of a transmitting device attached to an object to be located in a radiolocation scheme. The transmitter may have some Information Content that needs to be sent to the radiolocation network. For example the Information Content might be an intrusion detection sensor status. The information baseband is frequency translated to a frequency suitable to the radiolocation network application. A Carrier Frequency is applied to the I (In-phase) and Q (Quadrature-phase) Modulator that translates the Information Content to the desired radio frequency. The resultant signal has a finite bandwidth that is applied to, the input port of a common Mixer. The spreading signal is generated by the DSSS Source. The chip rate of the DSSS Source is determined by the Chip Clock Source. The frequency-spread signal is amplified and band limited by the Bandpass Amplifier and is launched into the air by the Antenna.

[0026] A network of receiving antennas receives the DSSS signal after the propagation of the signal through the air from the transmitter Antenna. Each receiving antenna gets the signal at times set by the respective distances between the Transmitter and individual Receivers.

[0027] A network of Radio Receivers is utilized to receive the combined signal. The combined signal at the receivers is created by the incident (line of sight) and the diffracted and reflected (non-line of sight) signals.

[0028] FIG. 2 illustrates a Receiver for a Radiolocation System for mitigating Correlator timing errors in the presence of Multipath Signals. FIG. 2 illustrates the significant components, including the Non-Linear Multipath Processing Device that improves the accuracy of locating a Transmitter in the presence of multipath signals. The combined signal arrives at a number of Antennas and is processed with the Receiver shown in FIG. 2. The Antenna transfers the combined signal to a Frequency Downconverter for convenience of processing the information. An Inphase and Quadrature Demodulator comprised of two diode mixers, a 90° Phase Shifter and a Local Oscillator. The signal is demodulated and decomposed into In-phase (I) and Quadrature-phase (Q) components. Each component, Q, and I is converted from an analog to a digital signal for time information extraction. A copy of the Transmitter's DSSS code word is sent to the Correlator for detection of the Transmitter's signal.

[0029] Quadrature Component Data (Q) data and Inphase (I) data is fed to the Correlators to extract the arrival time information. With only an incident signal the time measurement accuracy is limited to the processing circuitry and stability of local clocks. However the more typical case is when reflected signals add to the incident signal that the time measurement accuracy is limited by distortions created multiple copies of the incident wave delayed in time.

[0030] The problem becomes how to extract the time of arrival of the incident signal when corrupted with multiple path signals. The Non-Linear Multipath Processing Device estimates the composition of the incident and reflected signals and eliminates the measurement distortion to find the arrival time of the incident signal.

[0031] The corrected time information, from each radiolocation network receiver, is then sent to the Rectilinear Position Processor to determine the X, Y, and Z coordinates of the Transmitters location.

[0032] The Transmitter and Receiver clocks are not assumed to be phase-locked. That is, the transmitter carrier and receiver local oscillators may be out of phase. Each multipath signal is essentially the same as the original modulated carrier except for a phase shift, time shift, and amplitude change. We describe the combined signal received as follows:

[0033] m(t)=modulation as a function of time in seconds

[0034] ω=transmitter carrier frequency in rad/s

[0035] γ=phase of the receiver local oscillator in radians

[0036] β_k=amplitude of signal k, incident or multipath

[0037] θ_k=phase shift of signal k in radians

[0038] t_k=time delay of signal k in seconds  (0.1)

[0039] Then the received signal including the incident and all multipath signals is described by: $\begin{array}{cc}s\left(t\right)=\sum _k^{}{\beta }_{k}m\left(t-t_k\right)\cos \left(\omega \left(t-t_k\right)+{\theta }_k\right)& \left(0.2\right)\end{array}$

[0040] The time delays are assumed ordered by index. The propagation delay of the incident wave is t₀. The received signal is down-converted using the familiar trigonometric identity:

[0041] cos(x)cos(y)=½ cos(x+y)+½ cos(x−y)  (0.3)

[0042] The receiver local oscillator output is cos(ωt+γ) and cos(ωt+γπ/2) for an I/Q demodulator. The down-converted receiver signal is then: $\begin{array}{cc}I\left(t\right)=\sum _k^{}{\beta }_{k}m\left(t-t_k\right)\cos \left(\gamma +\omega t_k-{\theta }_k\right)Q\left(t\right)=\sum _k^{}{\beta }_{k}m\left(t-t_k\right)\sin \left(\gamma +\omega t_k-{\theta }_k\right)& \left(0.4\right)\end{array}$

[0043] We can represent these signals in complex form by associating the real part, as I and the imaginary part as Q. The coefficients of the time-shifted modulation are then:

[0044] α_k=β_k[cos(γ+ωt_k−θ_k)+j·sin(γ+ωt_k−θ_k)]  (0.5)

[0045] The complex valued baseband signal is then simply: $\begin{array}{cc}{s}_B\left(t\right)=\sum _{k=0}^{N_r}{\alpha }_{k}m\left(t-t_k\right),{\alpha }_k\mathrm{complex},N_r\mathrm{multipaths}& \left(0.6\right)\end{array}$

[0046] The number of significant multipath signals is N_r. A significant signal has a measurable effect on the receiver correlator output. What is needed is a technique to find N_r, {α_k} and {t_k}.

[0047] Correlation

[0048] Timing accuracy depends on the correlation properties of the reference signal. We assume that the reference signal is present in the received baseband signal, s_B(t), and attempt to correlate with it. Then, since correlation is linear, we have: $\begin{array}{cc}\mathrm{corr}\left(s_B,m\right)\left(t\right)=\sum _{k=0}^{N_r}{\alpha }_k\mathrm{corr}\left(m,m\right)\left(t-t_k\right)& \left(1.1\right)\end{array}$

[0049] Where the continuous correlation for real functions 𝑓 and m is defined as: $\begin{array}{cc}\mathrm{corr}\left(f,m\right)\left(t\right)={\int }_{-\infty }^{\infty }f\left(u+t\right)m\left(u\right)du& \left(1.2\right)\end{array}$

[0050] In practice, we approximate the correlation function with a discrete sum at the system Analog to Digital Conversion I/Q Sampling rate. The modulation signals, here m(·), is finite so that (1.2) is really a finite integral or sum. From equation (1.1), we see that the correlation of the received signal with the reference signal is a sum of time shifted copies of the autocorrelation function of the reference signal, corr(m,m)(t), with coefficients {α_k} and time shifts {t_k}. The information of interest is t₀, the time of arrival of the incident wave at the receiver.

[0051] The correlation of the received signal (without multipath but distorted by filtering in the transmitter and receiver) with the original modulating signal with be called autocorrelation. The correlation of the received signal including several possible multipath signals with the modulating signal will be referred to simply as correlation. That is, the sampled autocorrelation function is at the receiver where one of the reference signals originated as the modulation at the transmitter.

[0052] Splines

[0053] The autocorrelation function of the modulating signal is distorted due to filtering in the transmitter and receiver. This response is difficult to obtain analytically. In order to use non-linear filtering techniques later, we need to approximate the autocorrelation response at the receiver with a smooth curve that has a continuous derivative. To accomplish this, we use a least squares cubic spline. The data to fit is the sampled autocorrelation response at the receiver I/Q output.

[0054] A function 𝑓(x) is a cardinal spline of order n provided that:

[0055] 𝑓∈Cⁿ⁻²and 𝑓|_[k,k+2)∈π_n−1,k∈□  (2.1)

[0056] That is, 𝑓 has n−2 continuous derivatives and 𝑓 restricted to the interval [k,k+1) is a polynomial of degree at most n−1. This implies that the polynomial pieces match up smoothly at the integers, k∈□. The endpoints of these intervals are called knots. We will use n=4- cubic splines with 2 continuous derivatives and integer knots. A basis for the cardinal spline space of order n are the n^thorder cardinal B-splines, B_n(t). The cardinal B-splines can be defined explicitly as: $\begin{array}{cc}{B}_n\left(t\right)=\frac{1}{\left(n-1\right)!}\sum _{k=0}^n{\left(-1\right)}^k\left(\begin{array}{c}n\\ k\end{array}\right){\left(t-k\right)}_{+}^{n-1},\mathrm{where}{t}_{+}^{n-1}=\{\begin{array}{cc}{t}^{n-1}& t\ge 0\\ 0& \mathrm{otherwise}\end{array}& \left(2.2\right)\end{array}$

[0057] The cardinal B-splines may be defined recursively as:

[0058] B_n(t)=(B_n−1*B₁)(t)=∫₀¹B_n−1(u−t)du,B₁(t)=χ_[0,1)  (2.3)

[0059] Where χ_[0,1)is the characteristic function on the interval [0,1). The Cardinal B-splines have support [0,n] and are symmetric relative to their support as shown in FIG. 3.

[0060] FIG. 3 illustrates the Cubic Spline method of estimating a receiver correlation pulse response. Using this basis, a cardinal spline can be represented as a sum of centered and shifted cardinal B-splines as follows: $\begin{array}{cc}f\left(t\right)=\sum _{k=0}^N{c}_k{B}_n\left(t+n/2-k\right)& \left(2.4\right)\end{array}$

[0061] The derivative 𝑓′(t) will be needed later. Cubic splines have two continuous derivatives. The first derivative can be derived using the formula:

[0062] B′_n(t)=B_n−1(t)−B_n−1(t−1)  (2.5)

[0063] Equation (2.5) is substituted into (2.4) to yield a formula for generating the first derivative: $\begin{array}{cc}\begin{array}{c}{f}^{\prime }x=\sum _{k=0}^N{c}_k{B}_n^{\prime }\left(x+n/2-k\right)\\ =\sum _{k=0}^N{c}_k\left\{B_{n-1}\left(x+n/2-k\right)-B_{n-1}\left(x-1+n/2-k\right)\right\}\end{array}& \left(2.6\right)\end{array}$

[0064] The autocorrelation of the modulating signal (as measured at the receiver) is modeled using a least squares cubic spline as in (2.4) with n=4. A typical autocorrelation of a PN sequence is shown in FIG. 4.

[0065] FIG. 4 illustrates a Normal Autocorrelation Receiver output without received reflections. It is necessary to find the coefficients {c_k}, k=0, . . . ,N in equation (2.4) so that the spline function 𝑓(t) is approximately equal to the measured autocorrelation function corr(m,m)(t) at the sampled values. An assumption is made that there are 2^j·(N+1) total sample points with j≧1. Samples therefore occur at t=2^−j·i,i=0, . . . , N_swhere N_s=2^j·(N+1)−1. $\begin{array}{cc}\mathrm{corr}\left(m,m\right)\left(2^{-j}·i\right)\cong \sum _{k=0}^N{c}_k{B}_n\left(2^{-j}·i+\frac{n}{2}-k\right),i=0,\dots ,N_s& \left(2.7\right)\end{array}$

[0066] This is an over-determined system of linear equations. There are N_s+1 equations in the n+1 unknowns {c_k}, k=0, . . . , N. $a_{\mathrm{ik}}=B_n\left(2^{-j}·i+\frac{n}{2}-k\right),$

[0067] The (unweighted) least square solution is

[0068] ĉ=(A^TA)⁻¹A^Tz  (2.9)

[0069] The solution ĉ minimizes the residual sum of squares J(c) where

[0070] J(c)=(z−Ac)^T(z−Ac)  (2.10)

[0071] Once the coefficients {c_k} are computed, equations (2.4) and (2.6) are used to generate a table of ideal autocorrelation values and the first derivative of the ideal autocorrelation. The number of coefficients required in (2.4) is generally small as evidenced by FIG. 5 where only 5 coefficients are used with 32 sampled data points. FIG. 5 illustrates a Measured Received Signal and its Cubic Spline Approximation. The least squares cubic spline curve is superimposed on the measured autocorrelation curve to illustrate the accuracy of the approximation.

[0072] In practice, the sampled autocorrelation values are scaled to the resolution 2^−j. The tables for the ideal autocorrelation and derivative are generated using equations (2.4) and (2.6). The tables are then (implicitly) scaled back to the system sampling resolution.

[0073] The importance of the least squares cubic spline curve in generating the ideal autocorrelation and derivative tables is twofold:

[0074] The cubic spline has two continuous derivatives. The convergence of the iterative algorithm used later depends on the continuity of the first derivative.

[0075] Artificial local maximums or minimums caused by sampling the noisy receiver data at a high rate are removed. This prevents the iterative solution from being captured at a local max or min instead of converging toward the desired result.

[0076] Multipath Mitigation Algorithm

[0077] The proposed solution to the multipath problem is to find N_r, {α_k}, and {t_k} for k=0, . . . , N_r. Then the TOA is simply t₀. The algorithm used is a recursive iterated extended Kalman filter. The recursion is on the number of multipath signals, N_r=0,1,2, . . . . For each value of N_r, the non-linear equations involving {α_k}and {t_k}are solved by an iterated extended Kalman filter. Unique to this algorithm is the measurement noise model based on the process gain of a PN sequence and the termination criteria for N_rbased on a least square χ²statistic from the Kalman filter. In practice, a final value of N_r=2 is sufficient to resolve most multipath signals. The χ²statistic from the Kalman filter is used to generate a one-sided confidence interval (typically 95%) for the dimension N_rof the multipath model.

[0078] Correlation Noise Model

[0079] The measurement data in this problem is the sampled receiver correlator I/Q output. The correlator integrates over the time period of the PN sequence or other modulating signal. The autocorrelation of a single burst of a typical PN code has a two chip wide peak and multiple side peaks. The process gain of the PN code is defined as the ratio, of the amplitude in this peak to the amplitude of the next closest side peak. The gain is related to the length of the code. The measurement noise model used in the Kalman filter (see section 5.3) is as follows:

[0080] z=normalized measured Ior Q

[0081] σ_z=A_s·e^−k|z|+σ_Peak

[0082] A_s=amplitude of largest side peak

[0083] k=gain constant

[0084] σ_Peak=sample noise std. dev. at peak  (3.1)

[0085] The I/Q samples above a threshold determined by the receiver design and AGC State are first normalized. Each I/Q sample is a measurement z in the Kalman filter. The covariance of the measurement noise σ_z²is computed using equation (3.1). The gain k is chosen so that, σ_z≅σ_Peak, at the correlation peak. In practice, A_sand σ_Peakare a function of the receiver AGC state. FIG. 6 illustrates a Correlation Noise Peak Model. FIG. 6 illustrates the measurement noise model for a typical single correlation peak. The received signal has an SNR of 0 dB. The dotted line indicates the level A_s.

[0086] χ²Termination Statistic

[0087] The (weighted) least square problem is summarized in equation (3.2). Under the assumption that the measurement noise is normally distributed with a zero mean, it can be shown that the minimized least squares cost function J*=J({circumflex over (x)}) is χ²distributed with n_z−n_xdegrees of freedom. Note that cost function J(x) in equation (3.2) now includes the covariance of the measurement noise R in contrast to equation (2.10).

[0088] x—unknown state vector of dimension n_x

[0089] z—measurement vector of dimension n_z

[0090] H—linear measurement matrix of dimension n_z×n_x

[0091] v—noise vector of dimension n_z,v□N(0,R)

[0092] z=Hx +v

[0093] J(x)=(z−Hx)^TR⁻¹(z−Hx)  (3.2)

[0094] The Iterated Extended Kalman filter (described below) is used within the iterations without a priori information. The above statistic is therefore valid at the last iteration where H is the Jacobian of the non-linear measurement equation evaluated at the converged solution {circumflex over (x)}. If the statistic J* falls within the one-sided 95% confidence interval for a χ²distribution with n_z−n_xdegrees of freedom, the recursion on N_ris terminated. Note that a prerequisite is that the Iterated Extended Kalman filter converges at this step.

[0095] Kalman Filter

[0096] The output of the receiver correlator is the complex correlation of the superposition of the incident and all diffracted and reflected multipath signals with the modulating signal, typically a Pseudo Random Noise (PRN) code. The (complex-valued) correlation is modeled by equation (3.3). The correlation is the sum of time shifted sampled autocorrelation with complex amplitude coefficients {α_k} and continuous time shifts {t_k}, where k=0, . . . , N_r. The discrete sample times are {t_i}, i=1, . . . , N_s. $\begin{array}{cc}\mathrm{corr}\left(S_B,m\right)\left(t_i\right)=\sum _{k=0}^{N_r}{\alpha }_k\mathrm{corr}\left(m,m\right)\left(t_i-t_k\right),i=1,\dots ,N_s& \left(3.3\right)\end{array}$

[0097] The real and imaginary parts of the correlation are designated by I_iand Q₁as in equation (3.4). In practice, these are the normalized sampled receiver correlator outputs. $\begin{array}{cc}\begin{array}{c}{I}_i=\mathrm{Re}\left(\mathrm{corr}\left(S_{B,}m\right)\left(t_i\right)\right)\\ =\sum _{k=0}^{N_r}\mathrm{Re}\left({\alpha }_k\right)\mathrm{corr}\left(m,m\right)\left(t_i-t_k\right)\\ Q_i=\mathrm{Im}\left(\mathrm{corr}\left(S_{B,}m\right)\left(t_i\right)\right)\\ =\sum _{k=0}^{N_r}\mathrm{Im}\left({\alpha }_k\right)\mathrm{corr}\left(m,m\right)\left(t_i-t_k\right)\end{array}& \left(3.4\right)\end{array}$

[0098] To eliminate all complex arithmetic we use the definitions in equation (3.5):

[0099] α_Ik=Re(α_k)

[0100] α_Qk=Im(α_k)

[0101] x=(t₀, . . . ,t_N,α_I0, . . . ,α_IN,α_Q0, . . . ,α_QN)^T

[0102] h(x)=(I₁, . . . ,I_N,Q₁, . . . ,Q_N)^T

[0103] z_Ii=I_i+v_Ii,v_Ii=measurement noise

[0104] z_Qi=Q_i+v_Qi,v_Qi=measurement noise

[0105] z=(z_I1, . . . ,z_IN,z_Q1, . . . ,z_QN)^T  (3.5)

[0106] The non-linear measurement equation relative to the multipath model (3.3) is given by equation (3.6). Here the covariance of the measurement noise vector v is assumed diagonal with each diagonal entry computed as in equation (3.1) with z replaced by either I_ior Q_ifor i=1, . . . ,N_s.

[0107] z=h(x)+v, where v=(v_I1, . . . ,v_IN,v_Q1, . . . ,v_QN)^T

[0108] E[v_Ii]=0

[0109] E[v_Ii²]=(A_s·e^−k|I|+σ_peak)²

[0110] E[v_Qi]=0

[0111] E[v_Qi²]=(A_s·e^−k|Q|+σ_peak)²  (3.6)

[0112] The multipath effects can be mitigated by the solution of the non-linear measurement model presented by equations (3.3) through (3.6). The time shift t₀is then used as the time of arrival of the incident signal. The solution of this problem (for a fixed value of N_r) is the well-known iterated extended Kalman filter. Unique to the multipath problem are the measurement noise model as described in section (5.1) and the recursion on the number of multipath signals as described in section (5.2). The extended Kalman filter is a linearization of the non-linear measurement function (3.6) using a first order Taylor's series. The Jacobian of h(x) is required and computed in equation (3.7). $\begin{array}{cc}\begin{array}{c}H=\frac{\partial h}{\partial x}\\ =\left(\begin{array}{ccccccccc}\frac{\partial I_1}{\partial t_0}& \dots & \frac{\partial I_1}{\partial t_{N_r}}& \frac{\partial I_1}{\partial {\alpha }_{\mathrm{I0}}}& \dots & \frac{\partial I_1}{\partial {\alpha }_{{\mathrm{IN}}_r}}& 0& \dots & 0\\ & & & & ⋮& & & & \\ \frac{\partial I_{N_s}}{\partial t_0}& \dots & \frac{\partial I_{N_s}}{\partial t_{N_r}}& \frac{\partial I_{N_s}}{\partial {\alpha }_{\mathrm{I0}}}& \dots & \frac{\partial I_{N_s}}{\partial {\alpha }_{{\mathrm{IN}}_r}}& 0& \dots & 0\\ \frac{\partial Q_1}{\partial t_0}& \dots & \frac{\partial Q_1}{\partial t_{N_r}}& 0& \dots & 0& \frac{\partial Q_1}{\partial {\alpha }_{\mathrm{Q0}}}& \dots & \frac{\partial Q_1}{\partial {\alpha }_{{\mathrm{QN}}_r}}\\ & & & & ⋮& & & & \\ \frac{\partial Q_{N_s}}{\partial t_0}& \dots & \frac{\partial Q_{N_s}}{\partial t_{N_r}}& 0& \dots & 0& \frac{\partial Q_{N_s}}{\partial {\alpha }_{\mathrm{Q0}}}& \dots & \frac{\partial Q_{N_s}}{\partial {\alpha }_{{\mathrm{QN}}_r}}\end{array}\right)\end{array}& \left(3.7\right)\end{array}$

[0113] The first N_srows of the Jacobian relative to the I_isamples are computed easily in equation (3.8). $\begin{array}{cc}\begin{array}{c}\frac{\partial I_i}{\partial t_s}=\frac{\partial }{\partial t_s}\sum _{k=0}^{N_r}{\alpha }_{\mathrm{Ik}}\mathrm{corr}\left(m,m\right)\left(t_i-t_k\right)\\ ={\alpha }_{\mathrm{Is}}\frac{\partial }{\partial t_s}\mathrm{corr}\left(m,m\right)\left(t_i-t_s\right)\\ ={\alpha }_{\mathrm{Is}}\frac{d}{dt}\mathrm{corr}\left(m,m\right)\left(t\right){|}_{t=t_i-t_s}\\ \frac{\partial I_i}{\partial {\alpha }_{\mathrm{Is}}}=\mathrm{corr}\left(m,m\right)\left(t_i-t_s\right)\\ \frac{\partial I_i}{\partial {\alpha }_{\mathrm{Is}}}=0\end{array}& \left(3.8\right)\end{array}$

[0114] The last N_srows of the Jacobian relative to the Q_isamples are directly analogous and are computed in equation (3.9). $\begin{array}{cc}\frac{\partial Q_i}{\partial t_s}={\alpha }_{\mathrm{Qs}}\frac{d}{dt}\mathrm{corr}\left(m,m\right)\left(t\right){|}_{t=t_i-t_s}\frac{\partial Q_i}{\partial {\alpha }_{\mathrm{Is}}}=0\frac{\partial Q_i}{\partial {\alpha }_{\mathrm{Is}}}=\mathrm{corr}\left(m,m\right)\left(t_i-t_s\right)& \left(3.9\right)\end{array}$

[0115] The Taylor's series for h(x) expanded about x₀is given by equation (3.10). The Jacobian is evaluated using equations (3.7) through (3.9). $\begin{array}{cc}\begin{array}{c}h\left(x\right)\cong h\left(x_0\right)+{\left(\frac{\partial h}{\partial x}\right)}_{x_0}\left(x-x_0\right)\\ =h\left(x_0\right)+{\left(\frac{\partial h}{\partial x}\right)}_{x_0}x-{\left(\frac{\partial h}{\partial x}\right)}_{x_0}\mathrm{xo}\\ \mathrm{or}\\ h\left(x\right)-h\left(x_0\right)+{\left(\frac{\partial h}{\partial x}\right)}_{x_0}{x}_0={\left(\frac{\partial h}{\partial x}\right)}_{x_0}x\end{array}& \left(3.10\right)\end{array}$

[0116] In order to write equation (3.10) as a linearized Kalman filter equation, we rewrite the measurement equation as in (3.11) and substitute for h(x) in the (last line of) equation (3.10).

[0117] h(x)=z−v  (3.11)

[0118] The resulting equation, (3.12), is now a linear measurement equation with additive noise. Note that the actual measurement z is now replaced by the pseudo-measurement on the left side of equation (3.12). $\begin{array}{cc}z-h\left(x_0\right)+{\left(\frac{\partial h}{\partial x}\right)}_{x_0}={\left(\frac{\partial h}{\partial x}\right)}_{x_0}x+v& \left(3.12\right)\end{array}$

[0119] Note that all functions involved are based on the sampled autocorrelation and it's derivative. These values are interpolated from the tables computed in section (4). Linear interpolation is adequate. Interpolation is necessary since the values t_i−t_kare not necessarily in the table. $\begin{array}{cc}\frac{d}{dt}\mathrm{corr}\left(m,m\right)\left(t\right)& \left(3.13\right)\end{array}$

[0120] The multipath mitigation algorithm is now summarized. A threshold is set for the receiver correlator output. If |(I,Q)| exceeds this threshold, N_ssamples are taken and normalized to produce the measurement vector z as in equation (3.5). The number of multipath signals N_ris initialized to 0. The following three steps are performed until a solution is reached or the available processing is exceeded:

[0121] Step 1: Compute the Initial State Estimate, x₀. The times {t_k},k=0, . . . , N_rare equally spaced in the sample interval [t₁,t_N,]. The coefficients are initialized as α_1k=α_Qk={square root over (2/2)}, k=0, . . . , N_r.

[0122] Step 2: Compute the Kalman filter state estimate {circumflex over (x)} from equation (3.12), then compute the Relative Error |x₀−{circumflex over (x)}|/|x₀|. If the relative error is small, say 0.0005 or three digits; at this point convergence is achieved—go to step 3. Otherwise, set x₀={circumflex over (x)} and repeat step 2.

[0123] Step 3: Compute the χ²statistic J*=J({circumflex over (x)}) as in equation (3.2). If J* is in the 95% Confidence Interval for a χ²variable with n_z−n_xdegrees of freedom, then move onto step 4. Otherwise increment N_rand go to step 1.

[0124] Step 4: Accept t₀as the TOA.

[0125] If during any step the allotted time or resources is used up then the algorithm is terminated without giving an answer for the TOA. In this case the data associated with this receiver is discarded.

[0126] Simulation Examples

[0127] This algorithm was simulated with a 63-chip PN code with a chip rate of 20 MHz. The transmitter signal was BPSK modulated with a bandwidth of 20 MHz. The receiver bandwidth was set at 20 MHz and additive white Gaussian noise (AWGN) was added to the received signal to bring the SNR to 0 dB. Three random multipath signals were added with random phase, time delay and amplitude. The combined signal was demodulated and correlated with the PN code. The resulting correlator, Inphase (I) and Quadrature Phase (Q) outputs were then processed with the multipath mitigation algorithm. Several illustrative examples are now given.

[0128] FIGS. 7, 8 and 9 are from a signal with three distinct multipath components. The multipath mitigation algorithm decomposes the receiver correlator output into a superposition of time shifted sample autocorrelation as in equation (3.4). FIG. 7 illustrates the Inphase (I) components. The vertical line is at the TOA of the incident signal, which should coincide with the peak of the first I components. Each dotted peak represents an Inphase (I) component of the form α_1k·corr(m,m)(t−t_k).

[0129] FIG. 8 illustrates the equivalent Q components α_Qk·corr(m,m)(t−t_k). Note that the vertical line representing the Time Of Arrival of the incident signal is also at the peak of the first Quadrature Phase (Q) component since the Inphase (I) and Quadrature Phase (Q) components share the same time shift t_k.

[0130] FIG. 9 compares the original correlator I output (solid) with the mitigation algorithm model (dashed). Note that the fit is very good at the peaks. The noise model of section 5.1 emphasizes the peak data and de-emphasizes the noisy data near 0. The 0 dB SNR of the original signal is apparent away from the peaks.

[0131] FIG. 10 compares the correlator Q output (solid) with the mitigation algorithm Q model (dashed). Note that the multiple distorted peaks have been accurately modeled.

[0132] FIGS. 11, 12 and 13 are an example of shallow (close to the Transmitter being tracked) reflections or refraction partially canceling out the incident signal correlation peak. As before, the vertical line occurs at the incident signal correlation peak, which is, absent in both the Inphase (I) and Quadrature (Q) correlation. The model data (dashed) is a good fit to the correlator output (solid) even though only one multipath signal (N_r=1) is modeled. The first Inphase (I) and Quadrature Phase (Q) correlation peaks occur at the Time Of Arrival of the incident signal that occurs at the vertical line. Note that the first recognizable peak in the raw data would have produced a Time of Arrival error on the order of 100 nanoseconds (approximately 100 feet). The multipath mitigation algorithm as illustrated in FIG. 13 has salvaged this corrupted data.

[0133] FIG. 11 compares the correlator I output (solid) with the mitigation algorithm model (dashed). The vertical line is at the Time of Arrival of the incident wave, which does not produce a local peak. The model is a good fit to the data except near zero.

[0134] FIG. 12 compares the correlator Q output (solid) with the model Q data (dashed). The low amplitude of the raw data has degraded the fit. The vertical line occurs at the Time of Arrival of the incident signal. The Time of Arrival and correlation peak appears to be unrelated.

[0135] FIG. 13 presents both the Inphase (I) components (dotted) and the Quadrature Phase (Q) components (dashed). Note that the vertical line is at the correlation peak of the first I component. The corresponding Q component is barely visible. This is a normal result of the phase of the first I/Q component being near π. Note that the correct TOA is salvaged from the raw data by the mitigation algorithm.

[0136] The present invention has been described in relation to particular embodiments, which are intended in all respects to be illustrative rather than restrictive. Those skilled in the art will understand that the principles of the present invention may be applied to, and embodied in, various program modules for execution on differing types of computers and/or equipment, operating in differing types of networks, regardless of the application.

[0137] Alternate embodiments will become apparent to those skilled in the art to which the present invention pertains without departing from its spirit and scope. Accordingly, the scope of the present invention is described by the appended claims and supported by the foregoing description.