Blind despreading of civil GNSS signals for resilient PNT applications

Abstract

A method, system and apparatus are claimed for receiving, blindly despreading, and determining geo-observables, of true civil Global Navigation Satellite Systems (GNSS) navigation signals generated by any of the set of satellite vehicles and ground beacons, amongst false echoes and malicious GNSS signals from spoofers and repeaters; for identifying malicious GNSS signals, and preventing those signals from corrupting or capturing Pointing, Navigation, and Timing tracking operations; and for geolocating malicious GNSS signals. The invention also provides time-to-first-fix over much smaller time intervals than existing GNSS methods and can operate both in the presence of signals with much wider disparity in received power than existing techniques, and in the presence of arbitrary multipath. Further embodiments employing spatial/polarization diverse receivers that remove non-GNSS jammers received by the system, as well as targeted GNSS spoofers that can otherwise emulate GNSS signals received at victim receivers, are also claimed.

Claims

1. A blind, adaptive, despreading, and anti-spoofing reception system for blindly detecting, differentiating, and separating civil Global Navigation Satellite System (GNSS) signals generated by satellite vehicles (SVs), pseudolites, and nontargeted spoofers, without prior knowledge of or search over the GNSS ranging code for those signals, said system comprising: at least one antenna and at least one receiver, receptive to energy in the Global Positioning System (GPS) L1 band; identifying means for receiving, amplifying, and filtering at least one signal-in-space containing coarse acquisition (C/A) code from any of the set of GPS satellites, pseudolites, spoofers, and repeaters, with said identifying means connected to said antenna and receiver; at least one frequency converter that converts the at least one signal-in-space to a frequency converted signal having any of an intermediate and complex baseband frequency, connected to said means for receiving, amplifying, and filtering; at least one analog-to-digital convertor (ADC) that converts that frequency converted signal to any of a real data stream and complex data stream, connected to said frequency converter; digital signal processing hardware (DSP), connected to said analog-to-digital convertor; memory capacity connected to any of the set of the receiver, the identifying means, the at least one frequency converter, the at least one ADC, and DSP; and any of the set of software, firmware, and hardware any of installed on and connected with said DSP that performs blind adaptive processing operations on that at least one real or complex data stream.

2. A method to be used in a blind, adaptive, despreading, and anti-spoofing reception system as described in claim 1, said method further comprising: receiving with the identifying means at least one signal-in-space that is all of common to any of a set of globally-available and public signals in the GPS L1 band, contains coarse acquisition code, and purports to be from a SV within that system's field of view (FOV); downconverting that signal; noting both physical and logical aspects of what is received; implementing an adaptive despreading algorithm on the logical aspect to blindly detect, demodulate, and determine geo-observables of the C/A code; and, implementing comparative analysis of the physical aspects of the signal and physical and logical aspects of the coarse acquisition code, to determine whether the received signal is from a legitimate source and not a spoofer.

3. A method to be used in a blind, adaptive, despreading, and anti-spoofing reception system as in claim 2, wherein the step of implementing comparative analysis of the physical aspects of the C/A code to determine whether the received signal is from a legitimate source and not a spoofer, further comprises: examining the received signal by: demodulating that signal's navigation (NAV) signal element to produce a demodulated NAV result; comparing the demodulated NAV result against a standard which may be any of stored in and retrieved from at least one spoofer database to which the system has access, known from current ongoing processing, and obtained from internet sources; using that comparison to (a) determine a current, site-specific, Signal-to-Interference-and-Noise Ratio (SINR) norm for the signal, and (b) more finely determine any distortions in the frequency of acquisition, and thereby determine a current and local deduced SINR; comparing the deduced SINR and a value for the context-specific and time-specific SINR experienced by the unit; noting any discrepancies from the comparison; examining the received signal and using any detected preamble to resolve its Binary Phase-Shift Keying (BPSK) sign and subframe timing; first obtaining, and then using, the respective values for the SV ephemeris, timing, and Pseudo-Random Noise (PRN) from the NAV fields of the signal, comparing with a spoofer database these values to detect whether any of no, one, and many SV's are using the same PRN; evaluating the perceived time of acquisition (TOA) and frequency of acquisition (FOA) of the received signal in a Nyquist zone against the respective values held in the spoofer database of what should be present from the purported source; determining any of whether, when, and how these respective values match and fail to match; validating the received signal by using that determination; periodically updating its Positioning, Navigation, and Timing (PNT) tracker by using common geo-observables and SV data, comparing and updating any elements against and within the spoofer database, thereby correcting for tracker error degraded by new SVs; effecting a geolocation of any spoofer through analysis of any set of differences between the characteristics of the received signal and those respectively in the spoofer database; determining any set of its sample rate, sample phase, and Local Oscillator (LO) offset and updating with that any of its evaluative and self-positioning functions; and, excluding any signal evidencing any discrepancy for any of its geo-observable and ephemeris data comprising NAV, SINR, PRN, TOA, FOA, and PNT, from any further processing, thus reducing the computational load; for any validated received signal, effecting an active despreading of the received signal to process the logical, contained message therein.

4. A method to be used in a blind, adaptive, despreading, and anti-spoofing reception system as in claim 3, further comprising: classifying each validated received signal as a candidate signal; associating each candidate signal with its purported emitter; for each candidate signal, computing any set of the fine-TOA, FOA Nyquist zone, and DOA using an copy-aided parameter estimation implementation; determining whether, when, and how these respective values match and fail to match with known, experienced, and discoverable values for said candidate signal's purported emitter; for any non-validated signal, analyzing it to effect a geolocation of the emitter, now identified as a spoofer, based on the receiving system's known PNT data; and, for any validated signal, effect effecting an active despreading to process its message, engaging channelization implementations of any of multiple, alternative, comparative and cross-validating processing.

5. A method to be used in a blind, adaptive, despreading, and anti-spoofing reception system as in claim 4, wherein the step of downconverting that signal further comprises: converting the signal to a complex digital data stream with a complex sampling rate of M.sub.ADC samples/ms, where M.sub.ADC is an integer; and, passing the complex digital data stream through a 1:M.sub.ADC serial-to-parallel (S/P) converter, resulting in an M.sub.ADC×1 vector sequence with a 1 ksps (kilo samples per second) data rate.

6. A method to be used in a blind, adaptive, despreading, and anti-spoofing reception system as in claim 5, for multifeed systems where the signal is received over M.sub.feed spatially or polarization diverse feeds (one from each antenna), further comprising coherently converting each signal to complex baseband, and digitizing it at an M.sub.ADC ksps data rate, resulting in M.sub.feed data streams, each with an M.sub.ADC ksps data rate; then passing each of the data streams through a 1:M.sub.ADC S/P converter, and combining the output into an M.sub.DoF×1 vector with a 1 ksps data rate, where M.sub.DoF=M.sub.feedM.sub.ADC is the degrees of freedom (DoF) of the receiver.

7. A blind, adaptive, despreading, anti-spoofing system for radio transmission comprising: at least one antenna and at least one receiver, each receptive to energy in the GPS L1 band; connected to identifying means for receiving, amplifying, and filtering at least one signal-in-space containing coarse acquisition (C/A) code from any of the set of GPS satellites, pseudolites, spoofers, and repeaters; said identifying means further comprising a frequency converter that converts the at least one signal-in-space to a frequency converted signal having any of an intermediate and complex baseband frequency, connected to said means for receiving, amplifying, and filtering; said frequency converter further comprising at least one analog-to-digital convertor that converts that frequency converted signal to any of a real data stream and complex data stream; said analog-to-digital converter connected with digital signal processing hardware; memory capacity for storing any of the real data stream and complex data stream connected to both said digital signal processing hardware and said analog-to-digital convertor; and, at least one processor with any of installed software and connected firmware that performs blind adaptive processing operations on any of that real data stream and complex data stream.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0038] FIG. 1 is a drawing of the prior art, showing a model of the presumed normal condition wherein a receiver (1) is receiving GPS signals from a SV (3).

[0039] FIG. 2 is also a drawing of the prior art, showing the exemplary civil GNSS short-code signal generation model for Civil GNSS Transmission Operations.

[0040] FIG. 3 shows a tabular description of GNSS networks employing civil GNSS signals that are exploitable by the invention.

[0041] FIG. 4 is also a drawing of the prior art, showing an exemplary of the front-end reception processing for any single-feed receiver receiving GNSS signals; with 3(a) showing a direct-conversion receiver front-end that can be tasked to multiple GNSS bands; and 3(b) showing a front-end for a direct-conversion receiver that is dedicated to a specific GNSS band. Each can be used as a part of, and in, at least one embodiment of the invention, reducing any cost of both developing and deploying the invention.

[0042] FIG. 5 is also a drawing of the prior art, showing an exemplary of the front-end reception processing for any single-feed superheterodyne receiver receiving GNSS signals; with 4(a) showing a superheterodyne receiver front-end that can be tasked to multiple GNSS bands, and 4(b) showing a superheterodyne dedicated receiver front-end that is dedicated to a specific GNSS band. Each can be used as a part of, and in, at least one embodiment of the invention, reducing any cost of both developing and deploying the invention.

[0043] FIG. 6 is also a drawing of the prior art, showing the exemplary GNSS SV geo-observables (the metadata of any signal), including the position and transmit power, frequency, phase, and delay for any GNSS SV (1); the geo-observables for any GNSS receiver (1) including its position, orientation, and receive frequency, phase, and delay; and the signal's geo-observables, including its TOA, RIP (Received Incident Power), LOB (Line of Bearing), DOA, FOA, and CCP (Calibrated Carrier Phase).

[0044] FIG. 7 is a drawing of an exemplary scenario where a receiver (1) is receiving GPS signals from each of a SV (3) and a beacon (5) that is presumably on or near the ground. As before the receiver and SV have their geo-observables, as does the beacon (5) which also include its position and transmit power, frequency, phase, and delay; however, both the SV (3) and beacon (5) have for their respective signals their individual PRN.

[0045] FIG. 8 is a drawing of an exemplary, targeted spoofing scenario, where a receiver (1) receiving a signal from a SV (3) is also receiving a signal from a spoofer (7), that is attempting to deceive the receiver (1). The spoofer (7)

[0046] FIG. 9 is a drawing showing on the upper left the global orbital satellite system with the lines-of-transmission from a subset of 14 GNSS SV's (3a-3n) within the effective field of view of a target GNSS receiver at or near the surface of the globe, and on the bottom right depicting geometry local to that intended GNSS receiver (1), showing that the receiver is aboard a moving platform with altitude and velocities/accelerations consistent with a loitering high-altitude airborne unattended aerial system (UAS), and receiving transmissions from both the 14 GNSS SVs (3a-3n) and 10 additional terrestrial beacons (5a-5j), with altitudes and velocities/accelerations consistent with loitering small UAS's (sUAS's).

[0047] FIG. 10 is a collection of tabular details of summary signal TOA, FOA, and RIP parameters for the exemplary scenario of FIG. 6 measured in 1 millisecond increments for a 10 second interval.

[0048] FIG. 11 is a collection of tabular details of summary Azimuthal Signal and Elevation Signal LOB parameters, for the exemplary scenario of FIG. 6 measured in 1 millisecond increments for a 10 second interval.

[0049] FIG. 12 a collection of tabular details of summary Azimuthal Signal and Elevation Signal DOA parameters, for the exemplary scenario of FIG. 6 measured in 1 millisecond increments for a 10 second interval.

[0050] FIG. 13 shows the baseline single-feed receiver structure, describing the symbol-synchronous channelizer integral to the present invention.

[0051] FIG. 14 shows a specific instantiation of the symbol-synchronous channelizer, using Fast-Fourier Transform (FFT) operations. The general transform T_DoF and its FFT instantiation including a bin selection operation are used in ways unique and novel to the invention.

[0052] FIG. 15 shows the L1 GPS C/A Navigation Symbol Stream, and depicts attributes of that stream (decomposition into sum and difference components with unique features) that allow instantiation of blind despreading operations implemented in embodiments of the invention.

[0053] FIG. 16 shows an exemplary multifeed receiver with an optional beam-forming network (BFN) and an implementation of the symbol-synchronous channelization operations for the invention.

[0054] FIG. 17 shows an exemplary single-feed active despreader structure and implementation of the blind despreading for the invention.

[0055] FIG. 18 shows an exemplary multifeed active despreader structure and implementation of the blind despreading for the invention.

[0056] FIG. 19 shows a flowchart for the batch adaptation procedure including embodiments in which prior information, e.g., baseband navigation sequences or known transmitter frequencies-of-arrival and differences in FOA (DFOA), are, or can be made, available for the invention while processing the incoming transmission(s).

[0057] FIG. 20 shows a flowchart for the batch adaptation procedure for use in embodiments wherein the GPS C/A signal is being received by the invention while processing the incoming transmission(s) including embodiments in which prior information, e.g., baseband navigation sequences of known transmitter frequencies-of-arrival (FOA's) and differences in FOA (DFOA) are, or can be made, available for the invention while processing the incoming transmission(s).

[0058] FIG. 21 is a tabular summary of complexity of an exemplary instantiation of the blind despreader.

DETAILED DESCRIPTION OF THE INVENTION

Overview

[0059] The invention described herein is a blind, adaptive, despreading, anti-spoofing system for radio transmission comprising at least one antenna and at least one receiver, each receptive to energy in the GPS L1 band, connected to identifying means for receiving, amplifying, and filtering at least one signal-in-space containing coarse acquisition (C/A) code from any of the set of GPS satellites, pseudolites, spoofers, and repeaters; said identifying means further comprising a frequency converter that converts the at least one signal-in-space to a frequency converted signal having any of an intermediate and complex baseband frequency, connected to said means for receiving, amplifying, and filtering; with said frequency converter further comprising at least one analog-to-digital convertor that converts that frequency converted signal to any of a real data stream and complex data stream; and said analog-to-digital converter connected with digital signal processing hardware, memory capacity connected to said digital signal processing hardware for storing any of the real data stream and complex data stream and said analog-to-digital convertor; and at least one processor with any of installed software or connected firmware, that performs blind adaptive processing operations on any of that real data stream and complex data stream.

[0060] FIG. 10 through FIG. 12 summarize key TOA, FOA, RIP (FIG. 10), LOB (FIG. 11), and local DOA (FIG. 12) parameters obtained for the reception scenario shown in FIG. 9, measured in 1 millisecond (ms) increments over a 10 second (s) interval. Results are compiled separately for the GNSS SV's and the beacons, in order to contrast differences in statistics between the two transmitter types.

[0061] The SV and beacon signal TOA's adhere closely to a first-order drift model over 10 s intervals. However, the SV and beacon TOA's are distributed over distinctly different ranges, e.g., 68-86 ms for the SV's, and 0.21 to 1.01 ms for the sUAS beacons, consistent with 310 km ground horizon for a Reaper Unattended Aerial System (UAS, or high performance drone) flying at its operational altitude. Low-altitude platforms, e.g., small UAS's (sUAS's, e.g., commercial drones) and ground vehicular receivers will see beacon TOA's with even smaller ranges of value. In addition, the TOA drift rate is much lower for the beacons (0.16 μs/second median absolute drift rate) than for the SV's (1.3 μs/second median absolute drift rate), due to the much higher observed speed of the medium earth orbit (MEO) satellites.

[0062] The SV and beacon signal FDA's also adhere closely to a first-order drift model over 10 s intervals. However, the second-order difference is large enough to cause FOA shifts on the order of 1 Hz over a 10 s interval. Consequently, higher-order FOA models may be needed under some reception scenarios.

[0063] The SV and beacon signal RIP's adhere closely to a zero-order drift model over 10 s intervals, exhibiting a median first-order power drift of 0.00012 dB/s. However, the SV RIP's (−139 dBm median RIP) are much weaker than the beacon RIP's (−113 dBm median RIP), despite the 15 dB lower assumed transmit equivalent isotropic radiated power (EIRP, maximum directional energy) for the beacons.

[0064] The SV and beacon signal LOB's adhere to a zero-order azimuthal and elevation drift model over 10 s intervals, exhibiting a median absolute differential azimuth and elevation of 0.0076°/s and 0.0024°/s, respectively. However, the median absolute azimuth and elevation DOA drift rates jump to 0.31°/s and 0.012°/s, respectively. The azimuthal drift rate is nearly identical for all of the received signals, due to motion of the receiver, which completes a 360° circuit every 19 minutes in the scenario assumed here. The TOA and FOA drift will be substantive over reception intervals employed by invention. However, the effect of DOA drift will be minor if antennas that are isotropic in azimuth are employed by the receiver.

[0065] In an initial embodiment a single antenna receives a signal-in-space (SiS) and passes that through the identifying means that the system will use to frequency-convert the signal to complex baseband, and to convert the signal to a complex digital data stream with a complex sampling rate of M.sub.ADCf.sub.sym, where f.sub.sym is the symbol rate of the baseband data sequence being spread by the M.sub.chp-chip ranging code in FIG. 2, e.g., 1 ksps (kilosample per second) for any of the civil GNSS signals shown in FIG. 3, and M.sub.ADC is an integer. The integer M.sub.ADC is typically on the order of M.sub.chp, e.g., 1,023 for the GPS L1 C/A signal, in order to achieve the full processing gain of an adaptive despreader. However, it can be an arbitrary integer of convenience to the receiver. The signal is then passed through a 1:M.sub.ADC serial-to-parallel (S/P) converter, resulting in an M.sub.ADC×1 vector sequence with a data rate of f.sub.sym.

[0066] In alternative embodiments, the signal is received over M.sub.feed spatially or polarization diverse feeds (one feed for each antenna), coherently converted to complex baseband, and digitized at an M.sub.ADCf.sub.sym sampling rate, resulting in M.sub.feed data streams, each with an M.sub.ADCf.sub.sym data rate. In an optional further embodiment a beamforming network (BFN) is implemented between the antennas and the receivers. In this case, the antennas are coupled with M.sub.feed receivers through a BFN with an input port connected to each antenna, and with M.sub.feed output ports (generally not equal to the number of antennas, and typically less than the number of antennas). Each signal is then passed through a 1:M.sub.ADC serial-to-parallel (S/P) converter, and “stacked” into an M.sub.DoF×1 vector with a data rate of f.sub.sym, where M.sub.DoF=M.sub.feedM.sub.ADC is the degrees of freedom (DoF) of the receiver.

[0067] In theory, the resultant M.sub.DoF×1 signal sequence can be processed using fully-adaptive blind despreading methods that can achieve the maximum attainable signal-to-interference-and-noise ratio (SINR) of the despreader, allowing detection and demodulation of many more symbol sequences than conventional spread-spectrum receivers that use correlative or “matched filter” despreading methods, and successful demodulation at power differences that are much higher than the conventional jamming ratio of correlative despreaders. Moreover, blind despreading methods can acquire spread spectrum signals without a search over code family or timing offset, which can require acquisition times on the order of tens of minutes in cold-start scenarios.

[0068] Despite these advantages, however, the dimensionality of the linear combiner used in the adaptive despreader has generally entailed a computational processing cost that is too high for fully-adaptive blind despreading methods to be feasible or even possible. Gradient-based methods such as normalized least mean squares (NLMS) or normalized constant modulus algorithm (NCMA) with O(M.sub.DoF) complexity have been subject to extreme slow convergence in presence of strong jamming interference, due to the well-known wide eigenvalue spread condition that occurs in presence of strong interference. And rapidly-converging methods such as Affine Projections, block least-squares (BLS) and block least-square constant modulus algorithm (CMA; conjoined, BLS-CMA) have had prohibitively high O(M.sub.DoF.sup.2) complexity, and required adaptation over block sizes typically on the order of integer multiples of M.sub.DoF samples, e.g., multiple seconds, to be considered feasible for most uses or existing, non-specialized, commercially feasible hardware and firmware implementation(s).

[0069] However, as capabilities of modern digital signal processing equipment have continued to improve, fully adaptive blind despreading methods have become increasingly feasible. Motivated by these developments, the present invention and its implementation have been developed to realize the promise of these methods.

DETAILED DESCRIPTION OF THE DRAWINGS

[0070] FIG. 1, a drawing of the prior art, illustrates the presumed normal condition wherein at time T a receiver (1) is receiving GPS signals from a transmitter SV l (3). The receiver (1) has a position at time of transmission t represented by 3×1 vector p.sub.R (t) and an orientation at time of transmission t represented by 3×3 orientation matrix Ψ.sub.R(t); the SV l has a position at time of transmission t represented by 3×1 vector p.sub.T(t; l). From the point of view of the receiver (1), the observed position of SV l at time of transmission t is represented by 3×1 vector p.sub.TR(t; l)=p.sub.T (t; l)−p.sub.R(t).

[0071] FIG. 2, also a drawing of the prior art, shows a GNSS short-code signal generation model for Civil GNSS Transmission Operations, wherein for transmission (7) at link l, a symbol generator (112) takes an intended stream of navigation symbols b.sub.T (n.sub.NAV; l) (110) with sequence index n.sub.NAV and transmission rate f.sub.NAV, and generates a stream of transmit data symbols d.sub.T (n.sub.sym; t) with sequence index n.sub.sym and baseband symbol (transmission) rate f.sub.sym; while a code generator (113) takes a pseudorandom noise (PRN) code index k.sub.PRN(l) (111), taken from a known set of PRN codes, for that transmitter, and generates an M.sub.chp×1 PRN code vector c.sub.T(l)=[c(m.sub.chp; k.sub.PRN(l))].sub.m.sub.chp.sub.=0.sup.M.sup.chp.sup.−1. The stream of symbols d.sub.T(n.sub.sym; l) is combined with the spreading code at the mixer (multiplier element) (115), resulting in M.sub.chp×1 vector symbol sequence C.sub.T(l)d.sub.T (n.sub.sym; l), and fed to a M.sub.chp:1 parallel-to-serial (P/S) converter (116), resulting in a spread digital signal with rate M.sub.chpf.sub.sym. The spread digital signal is then passed to a digital-to-analog converter (DAC) (118), which pulse-amplitude modulates the digital signal at interpolation rate M.sub.chpf.sub.sym controlled by a transmitter clock, also operating at rate M.sub.chpf.sub.sym (117) (less adjustments to compensate for relativistic effects due to the SV's altitude and velocity), which controls both the DAC (118), and a local oscillator (LO) (119). The LO (119) generates a sinusoidal mixing signal at frequency f.sub.T(l), which combines with the DAC output signal in a transmitter (120) and transmitted through an antenna (121) to create a signal-in-space (SiS) comprising the complex baseband signal, frequency-shifted to transmit frequency f.sub.T(l) and with phase offset φ.sub.T(l) (122).

[0072] With the exception of the GLONASS (Russian) GNSS signal, the transmit frequency employed by each SV is identical, such that f.sub.T(l)≡f.sub.T. In addition, for all of the civil GNSS signals, f.sub.T(l) is an integer multiple of the baseband symbol rate, such that f.sub.T(l)=M.sub.T(l)f.sub.sym for some positive integer M.sub.T(l). However, neither assumption need hold, or even be known in some instantiations, for the invention to be implemented.

[0073] The interpolation function used in the DAC is also typically a rectangular pulse. However, this assumption need not hold, or even be known in some instantiations, for the invention to be implemented.

[0074] It should be noted that the processing steps shown in FIG. 2 is one conceptual means for generating civil GNSS signals, and ignore other SiS's that may be transmitted simultaneously with civil GNSS signals, for example, the GPS P, P(Y), M, and L1C SiS's. In practice, these signals are largely outside the band occupied by the civil GNSS signals, and moreover are received well below the noise floor. Therefore, their effect on the invention described herein is small. An important exception is military GNSS signals that may be generated by terrestrial pseudolites, as the non-civil components of those signals may be above the noise floor within the civil GNSS passband in some use scenarios.

[0075] The processing steps shown in FIG. 2 generate a particular class of spread spectrum signal, referred to herein as a modulation-on-symbol direct-sequence spread spectrum (MOS-DSSS) signal, described in more detail below, in which the spreading code used to spread each data symbol d.sub.T(n.sub.sym; l) is repeated within that data symbol, such that the rate M.sub.chpf.sub.sym spread data sequence generated by the P/S converter (116) can be expressed as c.sub.T(n.sub.chpmodM.sub.chp; l)d.sub.T(└n.sub.chp/M.sub.chp┘; l), where (•)modM and └•┘ denote the modulo-M operation and integer truncation operations, respectively. As shall be described below, this property induces inherent massive spectral redundancy to the civil GNSS format that is exploited by the invention described herein.

[0076] FIG. 3 provides a tabular description of GNSS networks employing civil GNSS signals that can be modeled MOS-DSSS signals. As this Figure shows, every GNSS network deployed to date possesses at least one signal type that can be modeled by FIG. 2. Moreover, all of these signals have a 1 millisecond (ms) ranging code period, and therefore a 1 ksps (kilosample per second) MOS-DSSS baseband symbol rate, albeit with widely varying chip rates and SiS bandwidths, allowing their demodulation using a common receiver structure. Lastly, as the last column of FIG. 3 shows, the symbol streams of each of these signals possess properties that can be used to detect and differentiate them from other signals in the environment, including military GNSS signals, jammers, and civil GNSS signals with the same or similar properties. For example, the GPS coarse acquisition (C/A) signal listed in the first row of FIG. 3 is constructed from a 50 bit/second (bps) BPSK (real) navigation sequence b.sub.T(n.sub.NAV; l), repeated 20 times per MOS-DSSS symbol to generate 1 ksps symbol sequence d.sub.T(n.sub.sym; l)=b.sub.T(└n.sub.sym/M.sub.NAV┘; l) where M.sub.NAV=20. Moreover, b.sub.T(n.sub.NAV; l) possesses internal fields that can be further exploited to aid the demodulation process, e.g., the TLM Word Preamble transmitted within every 300-bit Navigation subframe. FIG. 4, also a drawing of the prior art, shows exemplary direct-conversion front-end reception processing for any single-feed receiver receiving GNSS signals; with FIG. 4A showing a direct-conversion receiver front-end that is dedicated to a specific GNSS band, and FIG. 4B showing a direct-conversion receiver front-end that can be flexibly tasked to any GNSS band, consistent with a software defined radio (SDR) architecture. Each can be used as a part of, and in, at least one embodiment of the invention.

[0077] Both direct-conversion receivers shown in FIG. 4 comprise the following: [0078] An antenna (128) and low-noise amplifier (LNA) (131), which receives an SiS containing GNSS emissions within the FoV of the receiver. [0079] A local oscillator (LO) (133), which generates a complex sinusoid with frequency f.sub.R and phase φ.sub.R, comprising an in-phase (I) real sinusoid, and a quadrature (Q) sinusoid shifted by 90° in phase from the in-phase sinusoid; [0080] A complex mixer (134), which multiplies the LNA output signal by the conjugate of the complex sinusoid, thereby creating a complex signal with real and imaginary parts, referred to herein as the in-phase (I) and quadrature (Q) rails of that signal, which has been downconverted by frequency shift f.sub.R. [0081] A dual lowpass filter (LPF) (135), which filters each signal rail to remove unwanted signal energy from the complex downconverted signal, yielding analog scalar complex-baseband signal x.sub.R(t−τ.sub.R), centered at 0 MHz if f.sub.R=f.sub.T, where τ.sub.R is the delay through the receiver electronics. [0082] A dual analog-to-digital converter (ADC) (136), which samples each rail of x.sub.R(t−τ.sub.R) at rate f.sub.ADC=M.sub.ADCf.sub.sym, where f.sub.sym is the civil GNSS signal symbol rate shown in FIG. 2 and M.sub.ADC is a positive integer, thereby generating complex dual-ADC output signal x.sub.ADC(n.sub.ADC) with time index n.sub.ADC (137).

[0083] The receiver structure shown in FIG. 4A possesses an additional bandpass filter (BPF) (130), inserted between the antenna and LNA, and dedicated to a specific GNSS band, which removes unwanted adjacent-channel interference (ACI) from the received GNSS SiS, and depicts an LO that is locked to specific, preset value, consistent with a dedicated receiver that is optimized for a specific GNSS band. This structure is more typical of existing GNSS receivers, and provides enhanced signal quality at cost of receiver flexibility.

[0084] In contrast, the receiver structure shown in FIG. 4B is consistent with SDR receivers that can be tasked to any GNSS band under user control. This receiver structure can also be used to receive a segment of a GNSS band, e.g., a 1 MHz segment of the GPS L5 band, if f.sub.R is significantly offset from the GNSS band center f.sub.T. The ability to blindly detect, demodulate, and estimate geo-observables of any portion of any GNSS band, using the same reception, downconversion, and sampling hardware regardless of the actual bandwidth of the GNSS signals present in that band, is an important feature of the invention.

[0085] FIG. 5, also a drawing of the prior art, shows an exemplary of the front-end reception processing for any single-feed superheterodyne receiver receiving GNSS signals; with FIG. 5A showing a superheterodyne receiver front-end that is dedicated to a specific GNSS band, and FIG. 5B showing a superheterodyne receiver front-end that can be flexibly tasked to multiple GNSS bands. Each can be used as a part of, and in, at least one embodiment of the invention.

[0086] The heterodyne receivers shown in FIG. 5 shift the SiS received by an antenna (128) and amplified in an LNA (130) to IF frequency f.sub.IF, typically using a two-stage mixer as shown in this Figure, yielding analog scalar real-IF signal x.sub.R(t−τ.sub.R) centered on IF frequency f.sub.IF if f.sub.R=f.sub.T, with receive phase-shift φ.sub.R=φ.sub.R.sup.(1)+φ.sub.R.sup.(2) and delay τ.sub.R induced through the receiver electronics. The analog signal is then passed to an ADC (143) that samples x.sub.R(t−τ.sub.R) at rate M.sub.ADCf.sub.sym, where f.sub.sym is the civil GNSS signal symbol rate shown in FIG. 5 and M.sub.ADC is a positive integer, yielding real ADC output signal x.sub.ADC (n.sub.ADC) with time index n.sub.ADC (137). The receiver structure shown in FIG. 5A can be optimized for a specific GNSS band, while the receiver structure shown in FIG. 5B can be flexibly tasked to different GNSS bands, or to segments of GNSS bands, under software control.

[0087] As will be discussed below, the ADC output signal provided by either the direct-conversion or the superheterodyne receivers can be used by the invention, without any change to the invention design except control parameters used in channelization operations immediately after that ADC. The invention can also be used with many other receiver front-ends known to the art. In many cases, the invention is itself blind to substantive differences between those receiver front-ends.

[0088] FIG. 6 is also a drawing of the prior art, showing the exemplary observable parameters of signals transmitted from a GNSS SV (3) to a GNSS receiver (1) that can be used to geo-locate those signals (“geo-observables”) (6), as a function of parameters local to the GNSS SV (4) and the GNSS receiver (2). These geo-observables include the time-of-arrival (TOA) (after removal of delays induced in transmitter and receiver electronics), received incident power (RIP), line-of-bearing (LOB), direction-of-arrival (DOA) relative to the platform orientation, frequency-of-arrival (FOA) (after removal of known frequency offset between the transmitter and receiver center frequencies), and carrier phase (after calibration and removal of transmitter and receiver phase offsets), which are measurable and differentiable at the receiver (1).

[0089] FIG. 7 is a drawing of an exemplary scenario in which a receiver (1) is receiving GNSS signals from each of a SV (3) and a beacon (5) that is presumably on or near the ground. As before, the receiver can estimate measurable and differentiable observables for the beacon (13) and the SV (8). In addition, it is assumed that the beacon (5) and SV (8) are each employing spreading codes that are uniquely different, e.g., through use of PRN codes with different indices.

[0090] FIG. 8 is a drawing of an exemplary, targeted spoofing scenario, where a receiver (1) is receiving GNSS signals from each of a SV (3) and a spoofer (7), presumably on or near the ground, that has targeted that SV (3) and is attempting to deceive that receiver into replacing the SV's signal with its own signal. In particular, the spoofer is observing the targeted SV (3), and/or using prior knowledge of the targeted SV's position transmission parameters (4), and generating a spoofed signal with the same spreading code as the targeted SV, a transmit power that is strong enough to jam the targeted SV's received signal, and a delay and frequency shift that have been altered to closely match the TOA and FOA of the targeted SV, based on knowledge of the receiver parameters (14) that it has also obtained. Importantly, these parameters are not likely to include the DOA of the receiver (1), which is a function of its LOB and orientation (which it may not have accurate knowledge of) and the receive antennas location on the receiver platform (which it is likely to have no knowledge of). Nor are these parameters likely to include the calibrated carrier phase (CCP), which would require centimeter-level knowledge of the receiver position to spoof.

[0091] In some scenarios, the spoofer may be emulating the targeted SV (l′.sub.SV=l.sub.SV), e.g., in order to corrupt the navigation solution for the GNSS receiver, and may be attempting to closely its spoofed geo-observables to that SV. This is referred to as an “aligned” spoofing scenario. In other scenarios, the spoofer may be emulating a different SV than the targeted SV (l′.sub.SV≢l.sub.SV), in order to deceive the receiver into thinking that it is in another location entirely. This is referred to as a “nonaligned” spoofing scenario. In the latter case, multiple SV's in the field of view of the receiver must be spoofed. For practical reasons, all of these spoofed signals may be transmitted from the same location.

[0092] In some scenarios, the spoofer power level may be closely aligned with (but strong enough to jam) the power level of the targeted SV, in order to prevent the receiver from detecting the spoofing based on simple power-level differences between the SV and spoofer RIP. This is referred to as a “covert” spoofing scenario. In some scenarios, it may suffice to overwhelm the receiver with the spoofed SV signals, for example, if the spoofer is used in a “personal privacy spoofer (PPS)” to deceive a GNSS logger that is co-located with the PPS. Indeed, in such a case, it may not be feasible to closely align the RIP's of the spoofer and targeted SV signals.

[0093] FIG. 9 shows an exemplary reception scenario in which a loitering airborne GNSS receiver (1) is operating in the field of view of 14 GNSS SV's (3a-3n) and 10 loitering airborne beacons (5a-5j). The GNSS receiver dynamics are consistent with a Reaper unattended aerial system (UAS) flying at a 7,500 meter altitude, and moving along a 32.5 km diameter, 19 minute counter-clockwise circular track, consistent with a 89 meter/second (200 mps) constant airspeed and 1/20 gravity inward centripetal acceleration. The receive antenna is assumed to have isotropic gain, and is oriented in the direction of travel for the receiver, with a 2.7° inward roll consistent with the UAS acceleration. The GNSS SV's are assumed to be part of a 32-SV medium Earth orbit (MEO) GNSS network, organized into 8 orbital planes with 4 SV's per plane, placing between 7 and 16 SV's (median 12 SV's) within the FoV of the receiver over the course of a single day; 14 SV's are in the FoV of the receiver in this example. The terrestrial beacons' dynamics are consistent with small UAS's (sUAS's), flying at 12-to-120 meter altitudes, ±45 meter/second (±100 mph) airspeeds (negative airspeed for clockwise circular tracks), and 1/200 to 1/100 gravity inward centripetal acceleration. The SV's are each assumed to have a 15 dBW transmit EIRP, consistent with the GPS L1 C/A signals, and the beacons are each assumed to have a 0 dBW transmit power, consistent with commercially available low-power transmitters such as personal privacy devices (PPD's). This scenario provides an illustrative range of TOA's, FOA's, RIP's, LOB's, and DOA's for a reception scenario that is pertinent to the invention.

[0094] FIG. 10 provides a tabular summary of the TOA (FIG. 10A), FOA (FIG. 10B), and RIP (FIG. 10C) of the SV and beacon signals received under the reception scenario shown in FIG. 9, using statistics computed every millisecond over a 10 second interval. Results are compiled separately for the GNSS SV's and the beacons, in order to contrast differences in statistics between the two transmitter types and transmission platforms. Results are compiled separately for the GNSS SV's and the beacons, in order to contrast differences in statistics between the two transmitter types and transmission platforms. Parameters tabulated in this Figure include minimum, maximum, and median ranges; and results of linear regression analyses to determine closeness of fit of parameters to a zero-order time-variation model, in which parameters are constant over 10 seconds; a first-order time-variation model, in which parameters evolve linearly over 10 seconds); and a second-order time-variation model, in which parameters evolve at a quadratic rate over 10 seconds. Observations that can be made for this Figure are as follows: [0095] The SV and beacon signal TOA's adhere closely to a first-order TOA drift model over 10 s intervals, inducing a second-order differential TOA (DTOA) of <100 nanoseconds (ns) over 10 seconds. However, the SV and beacon TOA's are distributed over distinctly different ranges, e.g., 68-86 ms for the SV's, and 0.21-1.01 ms for the sUAS beacons, consistent with 310 km ground horizon for a Reaper UAS flying at 7.5 km. sUAS and ground vehicular receivers will see beacon TOA's with even smaller ranges of value. In addition, the first-order DTOA is much lower for the beacons (0.16 μs/second median absolute drift rate) than for the SV's (1.3 μs/second median absolute drift rate), due to the much higher observed speed of the MEO satellites. [0096] The SV and beacon signal FOA's also adhere closely to a first-order drift model over 10 s intervals, inducing a second-order differential FOA (DFOA) of <1 Hz over 10 seconds, in contrast to the >200 Hz difference in FOA between all of the platforms. However, the second-order DFOA is large enough to warrant higher-order FOA models under some reception scenarios. [0097] The SV and beacon signal RIP's adhere closely to a zero-order drift model over 10 s intervals, exhibiting a median first-order differential RIP (DRIP) of 0.00012 dB/s. However, the SV RIP's (−139 dBm median RIP) are much weaker than the beacon RIP's (−113 dBm median RIP), despite the 15 dB lower assumed transmit EIRP for the beacons.

[0098] FIG. 11 provides a tabular summary of the azimuthal (FIG. 11A) and elevation (FIG. 11B) LOB from the SV and beacon transmitters to the receiver under the reception scenario shown in FIG. 9, computed over the same time interval used in FIG. 10. The LOB is purely a function of the geometry between the transmitters and the receivers. Parameters tabulated in this Figure include minimum, maximum, and median ranges; and results of linear regression analyses to determine closeness of fit of parameters to zero-order, first-order, and second-order time variation models over the 10 s reception interval. As this Figure shows, the SV and beacon signal LOB's adhere to a zero-order azimuthal and elevation drift model over 10 s intervals, exhibiting a median absolute differential azimuth and elevation of 0.0076°/s and 0.0024°/s, respectively.

[0099] FIG. 12 provides a tabular summary of the azimuthal (FIG. 12A) and elevation (FIG. 12B) local DOA from the SV and beacon transmitters to the receiver under the reception scenario shown in FIG. 9, computed over the same time interval used in FIG. 10 and FIG. 11. The local DOA includes the effect of platform orientation on the observed directions of arrival from the transmitters, which varies from the LOB as the receiver platform moves along its flight path. Results are compiled separately for the GNSS SV's and the beacons, in order to contrast differences in statistics between the two transmitter types and transmission platforms. Parameters tabulated in this Figure include minimum, maximum, and median ranges; and results of linear regression analyses to determine closeness of fit of parameters to zero-order, first-order, and second-order time variation models over the 10 s reception interval. As this Figure shows, the receiver motion causes the local differential DOA to 0.3°/s in azimuth and 0.02° in elevation. In addition, the change in azimuthal differential DOA is nearly the same for both signal classes, indicating that the change in dynamics is almost completely due to motion of the receiver. However, in both cases, the DOA adheres strongly to a first-order model, i.e., the second-order differential change is low.

[0100] FIG. 13 shows a single-feed receiver structure used in one embodiment of the invention, describing the symbol-rate synchronous vector channelizer (150) integral to the present invention. A SiS is received at an antenna (128), and downconverted to complex baseband using a direct-conversion receiver (100) implementing operations such as those shown in FIG. 4, comprising reception, frequency downconversion by receive frequency f.sub.R to complex baseband representation, and sampling using a dual ADC (136) with sampling rate M.sub.ADCf.sub.sym, where M.sub.ADC is a positive integer. The sampled signal is then passed to a symbol-rate synchronous channelizer (150) comprising the following conceptual operations: [0101] A 1:M.sub.ADC serial-to-parallel conversion operation that transforms the rate M.sub.ADCf.sub.sym ADC output signal to an M.sub.ADC×1 vector signal x.sub.sym(n.sub.sym) with rate f.sub.sym, given by

[00001]$\begin{array}{cc}{x}_{\mathrm{sym}}\left(n_{\mathrm{sym}}\right)=\left(\begin{array}{c}{x}_R\left(T_{\mathrm{sym}}.Math.n_{\mathrm{sym}}-{\tau }_R\right)\\ \mathrm{.Math.}\\ x_R\left(T_{\mathrm{sym}}\left(n_{\mathrm{sym}}+\frac{m_{\mathrm{ADC}}-1}{M_{\mathrm{ADC}}}\right)-{\tau }_R\right)\end{array}\right).& \mathrm{Eqn}.Math..Math.\left(1\right)\end{array}$ [0102] An M.sub.ADC:M.sub.DOF cannelization operation (155), conceptually performed by multiplying x.sub.sym(n.sub.sym) by M.sub.DoF×M.sub.ADC channelization matrix T.sub.DoF (153) to form M.sub.DoF×1 channelizer output signal x.sub.DoF(n.sub.sym)=T.sub.DoFx.sub.sym(n.sub.sym) (159), where M.sub.DoF is the degrees of freedom (DoF's) of the channelizer output signal.

[0103] The data rate of x.sub.DoF(n.sub.sym) is nominally equal to the symbol rate of the GNSS signal of interest to the system (less offset due to nonideal error in the ADC clock), e.g., 1 ksps for the signals listed in FIG. 3. Hence this is referred to here as a symbol-rate synchronous vector channelizer. A variety of different channelization methods can be implemented using this general processing structure, including time and frequency domain channelization methods, or channelizations employing mixed time-frequency channelizations such as wavelet-based channelizers.

[0104] FIG. 14 shows a specific instantiation of the symbol-rate synchronous vector channelizer (150), using Fast-Fourier Transform (FFT) operations. The channelizer first applies a zero-padded FFT with M.sub.ADC input samples, M.sub.FFT≧M.sub.ADC output bins (152), and a (nominally rectangular) window {w.sub.FFT(m.sub.ADC)}.sub.m.sub.ADC.sub.=0.sup.M.sup.ADC.sup.−1 to each S/P output vector, where M.sub.FFT is conveniently chosen to minimize computational complexity and/or cost of the FFT operation. The channelizer then selects M.sub.DoF output bins {k.sub.bin(m.sub.DoF)}.sub.m.sub.DoF.sub.=0.sup.M.sup.DoF.sup.−1 for use in subsequent adaptive processing algorithms (154). The channelizer thereby implements channelizer matrix

[00002]$\begin{array}{cc}{T}_{\mathrm{DoF}}=\left[\exp \left(-j.Math..Math.2.Math.\pi .Math..Math.k_{\mathrm{bin}}\left(m_{\mathrm{DoF}}\right).Math.\frac{m_{\mathrm{ADC}}}{M_{\mathrm{FFT}}}\right).Math.w_{\mathrm{FFT}}\left(m_{\mathrm{ADC}}\right)\right].& \left(153\right)\end{array}$

[0105] In one embodiment, the output channelizer bins are preselected, e.g., to avoid known receiver impairments such as LO leakage, or to reduce complexity of the channelization operation, e.g., using sparse FFT methods. In other embodiments, the output channelizer bins are adaptively selected, e.g., to avoid dynamic narrowband co-channel interferers (NBCCI).

[0106] FIG. 14 is implemented for the real-ADC output of the superheterodyne receivers shown in FIG. 5, by passing the real ADC signal (137) into the same operations shown in this Figure, and choosing output bins centered on the intermediate frequency (IF) of the ADC output signal.

[0107] The approach extends easily to polyphase filters, given generally by x.sub.DoF(n.sub.sym)=

[00003]$\begin{array}{cc}\underset{m_{\mathrm{sym}}=0}{\overset{M_{\mathrm{sym}}}{.Math.}}.Math.T_{\mathrm{DoF}}\left(m_{\mathrm{sym}}\right).Math.x_{\mathrm{sym}}\left(n_{\mathrm{sym}}-m_{\mathrm{sym}}\right),e.g.,.Math.x_{\mathrm{DoF}}\left(n_{\mathrm{sym}}\right)=T_{\mathrm{DoF}}.Math.\underset{m_{\mathrm{sym}}=0}{\overset{M_{\mathrm{sym}}}{.Math.}}.Math.\left(w_{\mathrm{sym}}\left(m_{\mathrm{sym}}\right)\odot x_{\mathrm{sym}}\left(n_{\mathrm{sym}}-m_{\mathrm{sym}}\right)\right),& \mathrm{Eqn}.Math..Math.\left(2\right)\end{array}$

where {w.sub.sym(m.sub.sym)}.sub.m.sub.sym.sub.=0.sup.M.sup.sym.sup.−1 is an order-M.sub.sym polyphase filter, “⊙” denotes the element-wise multiplication operation, and

[00004]$T_{\mathrm{DoF}}=\left[\exp \left(-j.Math..Math.2.Math.\pi .Math..Math.k_{\mathrm{bin}}\left(m_{\mathrm{DoF}}\right).Math.\frac{m_{\mathrm{ADC}}}{M_{\mathrm{ADC}}}\right)\right]$

is the M.sub.DoF×M.sub.ADC matrix that implements an unpadded, unwindowed FFT with M.sub.ADC input samples and M.sub.DoF output bins. This channelizer is especially useful in environments subject to very strong NBCCI. However, it affects dispersion added to the channelizer output signal.

[0108] The general transform T.sub.DoF and its FFT instantiation including the bin selection operation, are used in ways unique and novel to the invention.

[0109] FIG. 15 shows the L1 GPS C/A Navigation Symbol Stream, d.sub.T(n.sub.sym; l) (160), and depicts attributes of that stream (decomposition into sum and difference components {tilde over (d)}.sub.T(n.sub.sym; l).sub.+ (161) and {tilde over (d)}.sub.T(n.sub.sym; l).sub.− (162), respectively, with unique features) that allow instantiation of blind despreading operations implemented in embodiments of the invention. As this Figure shows, {tilde over (d)}.sub.T(n.sub.sym; l).sub.+ (161) only has zero value at transitions between navigation bits, whereas {tilde over (d)}.sub.T(n.sub.sym; l).sub.− (162) only has nonzero value at those transitions. The signal components also possess distinctively different features, which can be used to both detect and differentiate between those signal components in subsequent adaptive processing operations.

[0110] FIG. 16 shows an exemplary multifeed extension of the single-feed receiver, in which the GNSS signals are received over multiple spatial or polarization diverse antennas. The multifeed receiver front-end and symbol-rate synchronous channelization operations start when GNSS signals are received on multiple antennas (90a, 90b, 90c) using an array of spatial and/or polarization diverse antennas. In one embodiment, there are M.sub.feed antennas, and the signals received by the array are passed directly to a bank of direct-conversion receivers (101) where they are coherently down-converted to M.sub.feed complex-baseband signals, represented here as M.sub.feed×1 vector signal x.sub.R(t−τ.sub.R)=[x.sub.R(t−τ.sub.R; m.sub.feed)].sub.m.sub.feed.sub.=1.sup.M.sup.feed. In another embodiment, the antenna output signals (90a, 90b, 90c) are passed to a beamforming network (BFN) (94), which converts those signals to M.sub.feed signals that are in turn passed to the receiver bank (101). The receiver output signal vector is then sampled by a bank of M.sub.feed, dual ADC's (143), each with sampling rate M.sub.ADCf.sub.sym, provided by a common clock (117), where M.sub.ADC is a positive integer and f.sub.sym is the baseband GNSS signal symbol rate, e.g., 1 ksps for all of the civil GNSS signals shown in FIG. 3. The M.sub.feed ADC output signals are then passed to a bank of symbol-rate synchronous vector channelizers (160), each synchronous vector channelizer implemented using a 1:M.sub.ADC S/P operation and a M.sub.chn×M.sub.ADC channelization matrix T.sub.chn, in one embodiment, which converts the signals into M.sub.feed M.sub.chn×1 signals {x.sub.chn(n.sub.sym; m.sub.feed)}.sub.m.sub.feed.sub.=1.sup.M.sup.feed. The channelizer outputs are then combined into a M.sub.DoF×1 complex vector representation x.sub.DoF(n.sub.sym) (159), for example, by stacking the channelized signals as

[00005]$\begin{array}{cc}{x}_{\mathrm{DoF}}\left(n_{\mathrm{sym}}\right)=\left(\begin{array}{c}{x}_{\mathrm{chn}}\left(n_{\mathrm{sym}};1\right)\\ \mathrm{.Math.}\\ x_{\mathrm{chn}}\left(n_{\mathrm{sym}};M_{\mathrm{feed}}\right)\end{array}\right),& \mathrm{Eqn}.Math..Math.\left(3\right)\end{array}$

where M.sub.DoF=M.sub.feedM.sub.chn is the total number of receiver DoF's.

[0111] The multifeed receiver structure shown in FIG. 16 is one of many structures that the invention can be used with. For example, the bank of direct-conversion receivers (101) and dual-ADC's (143) can be replaced by a bank of superheterodyne receivers and real ADC's, and the channelizer bank (160) can be altered to convert the subsequent real-IF ADC output vector to the same complex channelizer output signal Note that the structure of the symbol-rate channelized multifeed signal x.sub.DoF(n.sub.sym) (159) shown in FIG. 16 is identical to the structure of the single-feed x.sub.DoF(n.sub.sym) (159) shown in FIG. 13. In fact, as shall be shown below, the blindly signal detection, demodulation, and geo-observable estimation methods used in the invention can be applied to either a single-feed or a multifeed receiver with no change in structure, albeit with differences in implementation complexity if M.sub.DoF varies as a function of M.sub.feed, and with differences in performance, e.g., if the additional spatial/polarization of the received signals allows the processor to better separate GNSS signals or remove co-channel interference (CCI), e.g., jammers.

[0112] To overcome this complexity issue, in one embodiment, the baseline channelizer method is adjusted to provide the same number of DoF's for the multifeed receiver as the single-feed receiver, for example, by reducing M.sub.chn or the ADC sampling rate by 1/M.sub.feed, such that M.sub.DoF remains constant as M.sub.feed is increased. This is easily accomplished using the thinned-FFT channelizer shown in FIG. 14, by simply reducing the number of output bins used by the channelizer. If the channelizer bins remain densely spaced, this approach has the added benefit of reducing the effect of TOA blur on the output signal. Alternately, if the channelizer bins are sparsely spaced, this approach can reduce complexity of the channelization operation itself. However, it should be noted that this processor sacrifices despreading gain that could improve the SINR of the output symbol stream.

[0113] Importantly, while the TOA and FOA of a GNSS transmitter can be easily spoofed in the targeted spoofing scenario shown in FIG. 8, the DOA (and, to a lesser degree, the polarization) of that transmitter cannot be easily spoofed. In addition, the multifeed receiver can null any CCI impinging on the array, if the array has sufficient degrees of freedom to separate that CCI from the GNSS signals. This includes both targeted spoofers and wideband jammers.

[0114] FIG. 17 shows an exemplary single-feed despreading structure. The GNSS SiS is received by an antenna (128) and passed to a receiver front-end, comprising at least a receiver and an ADC converter (170) with sampling rate M.sub.ADCf.sub.sym, where M.sub.ADC is a positive integer and f.sub.sym is the baseband symbol rate of the civil GNSS signals, followed by a symbol-rate synchronous vector channelizer (171), which downconverts, digitizes, and channelizes the GNSS SiS to provide an M.sub.DoF×1 vector signal x.sub.DoF(n.sub.sym) with sample rate f.sub.sym equal to the baseband symbol rate of the GNSS signals. The symbol-rate synchronous channelizer (172) passes the vector signal on to: [0115] A blind despreader (180), which detects, despreads, and demodulates the received GNSS signal without explicitly using or searching over ranging codes, using only the MOS-DSSS properties of the civil GNSS signals. Among other advantages, this can operate under potentially heavy CCI. This stage also estimates at least the fractional fine FOA (less a factor-of-1 kHz ambiguity) and DFOA of each detected signal, and the coarse (1 ms granularity) TOA of each detected signal. In some embodiments, this stage also estimates the MIMO channel signature, described below, for each detected signal. [0116] If desired, the blind despreader (180) passes these estimates to a copy-aided analyzer (184), which determines remaining geo-observables needed for a full PNT solution, using copy-aided parameter estimation methods that exploit results of the blind despreading stage and detailed structure of the MIMO link matrices, and ranging code information held in memory (174). Although this stage does require search over code timing to determine fine TOA, the search is simplified by geo-observables already estimated in the blind despreader, and does not require a search over code PRN if the PRN index is provided in the demodulated navigation signal stream. The copy-aided estimator also outperforms correlative methods, especially in presence of strong CCI and GNSS signals, e.g. co-channel beacons, due to inherent robustness of the copy-aided approach to CCI and network interference. Lastly, given calibration of the transmitter and receiver phase, the copy-aided analyzer also provides an estimation of the carrier phase suitable for high-precision geolocation.

[0117] The copy-aided analyzer then passes any of the set of demodulated SV ID's, ephemeres, observed TOA, FOA, DFOA, CCP, RIP—the geo-observables—for each detected signal to a pointing, navigation and timing (PNT) computation element (190), implemented using methods well known to those skilled in the art. In some embodiments, the PNT computation element (190) further sends timing and clock offset adjustment parameters back to the receiver front-end (170), in order to remove any FOA or TOA offset due to error between the receiver and GNSS transmitter clocks.

[0118] FIG. 18 shows an exemplary multifeed extension of the despreader structure shown in FIG. 17. A GNSS SiS is received by an array of spatially and/or polarization diverse antennas (129a-129N), and passed to a receiver bank (172), comprising M.sub.feed receivers and M.sub.feed ADC converters, each with sampling rates M.sub.ADCf.sub.sym, where M.sub.ADC is a positive integer and f.sub.sym is the baseband symbol rate of the civil GNSS signals, followed by a bank of symbol-rate synchronous vector channelizers (173) with output channelizer dimensionality M.sub.chn, which coherently downconverts, digitizes, and channelizes the GNSS SiS to provide an M.sub.DoF×1 vector signal x.sub.DoF(n.sub.sym) with sample rate f.sub.sym, where M.sub.DoF=M.sub.feedM.sub.chn.

[0119] The channelizer output signal is then passed to a blind despreader (180), which implements the same operations performed in FIG. 17, and a copy-aided analyzer, which determines remaining geo-observables needed for a full PNT solution, using copy-aided parameter estimation methods that exploit results of the blind despreading stage and detailed structure of the multiple-input, multiple-output (MIMO) link matrices, and ranging code information and array manifold information held in memory (174). Using array manifold information, the analyzer is able to further analyze this information to determine the DOA of the detected signals.

[0120] The copy-aided analyzer then passes any of the set of demodulated SV ID's, ephemeres, and geo-observables (including DOA) for each detected signal to a pointing, navigation and timing (PNT) computation element (190), implemented using methods well known to those skilled in the art. In some embodiments, the PNT computation element (190) further sends timing and clock offset adjustment parameters back to the receiver front-end (170), in order to remove any FOA or TOA offset due to error between the receiver and GNSS transmitter clocks.

[0121] Among other advantages, if coupled to a spatial/polarization diverse antenna array as shown in FIG. 18, the blind despreader provides additional inherent excision of non-GNSS CCI impinging on the array, including CCI from up to M.sub.feed−1 noiselike jammers. The blind despreader also provides inherent excision of targeted spoofers impinging on the array, even if those spoofers employ the same PRN code as the targeted GNSS signal and are closely aligned in FOA and TOA with that targeted GNSS signal, and even if their RIP is close to the RIP of the targeted GNSS signal (covert aligned spoofing scenario) or much stronger than the target GNSS signal (noncovert aligned spoofing scenario), due to the differing DOA of the spoofer and targeted GNSS signal. This advantage is not provided through any additional processing performed in the blind despreader, but simply through inherent properties of the transmission channel under this reception scenario.

[0122] The copy-aided analyzer shown in FIG. 18 can also be implemented using joint TOA-DOA estimators that are substantively similar to the TOA-estimators employed in the single-feed system shown in FIG. 17, or they can be implemented using methods the explicitly exploit channel structure introduced by the multifeed front-end. In some use scenarios where it suffices to simply identify SV's and ground beacons, the need to search over TOA and DOA can be eliminated entirely (using methods discussed in a further alternative and novel embodiment—not shown).

[0123] The despreading structures shown in FIG. 17 and FIG. 18 can be further simplified if prior information such as trial signal FOA's or GNSS navigation streams are made available to the receiver (using methods discussed in a further alternative and novel embodiment, e.g., through web-based ancillary data capture—not shown).

[0124] FIG. 19 shows a flowchart for batch processor, applicable to the single-feed or multi-feed despreader structures shown in FIG. 17 and FIG. 18, in which channelizer output data is collected over N.sub.sym output samples and processed on a stand-alone basis to detect and extract the GNSS data contained within that sample set. This flow diagram encompasses embodiments in which prior information, e.g., GNSS navigation sequences or known transmitter frequencies-of-arrival (FOA) and differences in FOA (DFOA), are, or can be made, available for the invention while processing the incoming transmission(s).

[0125] The first stage (blind despreading component) of this process comprises the following key processing operations: [0126] A channelized data whitening operation (220), described in more detail below, in which the data block is separated into a whitened data component and a set of autocorrelation statistics using highly regular and well-known linear algebraic operations, e.g., QR decompositions, implemented as a means for reducing complexity of subsequent processing steps. [0127] An emitter detection substage (221), which detects received data modes due to presence of GNSS signals in the environment, using known properties of the underlying symbol sequence, and provides preliminary estimates of the fine fractional FOA's and DFOA's of those modes. In presence of TOA blur due to nonzero DTOA, each GNSS signal can generate multiple modes due to TOA drift over the batch reception interval. [0128] A FOA refinement substage (223), in exemplary embodiments performed as part of the emitter detection stage, which further refines fine FOA's and DFOA's of modes associated with the GNSS signals. [0129] A mode association substage (225), which associates those modes with specific GNSS transmissions based on commonality of the FDA's and DFOA's for those modes. [0130] A despreader refinement substage (227), which refines the despreading solution and fine FOA and DFOA estimates of the detected GNSS signals based on additional properties of the underlying symbol sequence. [0131] A navigation (NAV) stream demodulation substage (228), which identifies the navigation stream timing phase (coarse TOA) within each GNSS symbol stream, and extracts the navigation symbols from the symbol stream. [0132] Optional despreader and FOA refinement substages (229) that further refine despreading solutions, fine FOA, DFOA, and coarse TOA of the detected and demodulated GNSS signals, e.g., using decision feedback operations. [0133] An optional MIMO link signature estimation substage (231), which blindly estimates the MIMO link signatures of the detected GNSS signals, for use in subsequent copy-aided analysis stages.

[0134] If needed, the spatial signatures and NAV data are then passed to a copy-aided DOA estimation stage (230), implemented if the invention is implemented using a multifeed receiver such as that shown in FIG. 18, which uses array manifold information stored in memory (230) to estimate the direction of arrival of the detected emitters; and to a copy-aided TOA estimation stage (241), which uses ranging code information stored in memory (232) to estimate the fine (sub-symbol) TOA, FOA Nyquist zone (multiple of 1 kHz FOA ambiguity), and carrier phase of the detected GNSS symbols. In some embodiments, the DOA and TOA estimators can be decoupled, or simplified estimators that first determine coarse DOA estimates of the detected signals can be employed to quickly identify terrestrial emitters, including beacons and spoofers, based on elevation of the emitters (both emitter classes), or commonality of DOAs (e.g., spoofers emanating from a single source).

[0135] Processing stages (221), (223) and (225) are typically only performed in “cold start” scenarios (218) where no information about the GNSS signals in the received data band are available to the receiver. Example use scenarios include those where the receiver has just been turned on; where the receiver has been idle for an extended period (co-called “warm start” scenarios); or where the receiver has tuned to a new band to assess availability of GNSS signals within that band. These processing stages can be eliminated or simplified if prior information can be made available to the receiver, e.g., known FOA's or NAV data, through a memory database (222).

[0136] These approaches are referred to here as fully-blind and partially-blind despreading strategies. In both approaches, the channel signature matrix, and therefore the ranging codes underlying that signature matrix, are assumed to be unknown (e.g., “Type C3” blind reception strategies as defined in the prior art). This assumption allows the use of broad ranges of methods that can operate without prior knowledge of the GNSS ranging code, or elements of the channel or receiver that can modulate that code, including FOA, TOA, or multipath induced by the channel; frequency, sampling, and timing offsets induced by the receiver front-end; and potentially severe filtering induced by the receiver. This assumption also allows the use of methods that can approach or achieve the maximum-attainable SINR (max-SINR) of the despreader, using linear algebraic methods that do not require complex prior information about the GNSS signals or interferers in the field of view of the receiver.

[0137] Neither ranging codes (232) nor (for multifeed receivers) array manifold data (230) is required until the copy-aided TOA (241) and DOA (240) analysis stages. In particular, the navigation data can be demodulated without ever searching over that data. This is a marked departure from both conventional correlative or despreading methods, which require identification and search over the ranging code before the navigation signal can be demodulated, and “codeless despreading” algorithms that cannot demodulate that data.

[0138] FIG. 20 shows a flowchart for a batch adaptation procedure for use in embodiments wherein the GPS C/A signal is being received by the invention, including embodiments in which prior information, e.g., baseband navigation sequences of known transmitter frequencies-of-arrival (FOA's) and differences in FOA (DFOA) are, or can be made, available for the invention while processing the incoming transmission(s).

[0139] The first stage (blind despreading component) of this process comprises the following key processing operations: [0140] A channelized data whitening operation (220), in which the data block is separated into a whitened data component and a set of autocorrelation statistics using highly regular and well-known linear algebraic operations, e.g., QR decompositions, implemented as a means for reducing complexity of subsequent processing steps. [0141] An emitter detection operation (251), which detects GNSS emitters in the environment using a phase-SCORE algorithm such as those described in the prior art, that exploits the high single-lag correlation of the 1 ksps baseband GNSS symbol streams. [0142] A FOA vector refinement operation (253), which refines the FOA vectors provided by each detected phase-SCORE mode, using a generalized frequency doubler spectrum that exploits the BPSK structure of the GNSS symbol streams, and prunes the detected modes based on statistics generated during the FOA vector optimization procedure. [0143] A mode association operation (255), which associates confirmed phase-SCORE modes with specific emitters based on commonality of the FOA vectors for modes induced by each emitter. [0144] A further despreader refinement procedure (256), which refines FOA vectors and generates refined 1 ksps FOA-compensated symbol data, using a generalized conjugate-SCORE algorithm. [0145] A navigation stream demodulation algorithm (257), which estimates the coarse (1 ms granularity) TOA of the despread 1 ksps symbol data and extracts the 50 bps NAV sequence from that data. [0146] If sufficient data is available, analyze the demodulated navigation stream to detect the telemetry (TLM) preamble transmitted every 6 seconds within the NAV stream, and use that preamble to resolve ambiguities induced by the blind adaptation procedures (±1 BPSK ambiguity, modulo-20 TOA symbol ambiguity), and to extract additional navigation information needed to complete a navigation solution for the receiver, e.g., the handover word (HOW) message, and the SV ephemeres (258). [0147] Further refine parameters based on the demodulated NAV message, e.g., using maximum-likelihood (MI) estimation (MLE) strategies (259). [0148] Estimate MIMO link signatures (231) for use in subsequent copy aided DOA (240) and copy-aided TOA (241) operations, if needed for use scenarios that the invention is being applied to.

[0149] As in the more general flow diagram shown in FIG. 19, steps (251), (253), and (255) are only needed in cold start scenarios (218) where no information about the GPS signals (except perhaps the GPS almanac) are available to the receiver. If additional material can be made available, then these steps can be dispensed with. In particular, if prior FOA's and DFOA's can be made available to the receiver (250), then processing can start with the generalized conjugate SCORE procedure (256). If the 50 bps NAV stream can be made available to the receiver (260), without or without prior FOA's, processing can start at the maximum-likelihood (ML) estimator (259).

[0150] FIG. 21 shows tabulated complexity of the blind despreader, for a fully-blind instantiation applied to the GPS C/A Signal. The Figure assumes reception of 24 civil GNSS signals, e.g., 12 signals from GNSS SV's and 12 signals from terrestrial beacons or spoofers, using a four-feed receiver (M.sub.feed=4), the FFT-based symbol-rate synchronous vector channelizer shown in FIG. 14, with 16, 32, 64, 128, or 256 output FFT bins/feed, and the blind despreader shown in FIG. 20, with 64, 128, 256, 512, or 1,024 degrees of freedom. Complexity results are summarized at the end of the phase-SCORE detection operations, which includes complexity needed to channelized and whiten the ADC output data, and to implement all of the phase SCORE detection operations using that whitened data; and at the end of the conjugate SCORE refinement operations. This Figure does not directly address the performance of the despreader in any of these processing scenarios; however, all of these scenarios should be able to separate all 24 signals, and all of these scenarios should in particular be able to detect and demodulate the terrestrial signals if they are emanating from beacons with reasonable (e.g., 30 dBm) transmit EIRP's. Moreover the scenarios with 256 DoF's should be able to detect and demodulate any signals received with RIP's on the order of −140 dBm or higher, commensurate with SV and beacon/spoofer RIP's under nominal reception conditions.

[0151] Observations that can be made from this Figure are as follows: [0152] The channelization and whitening operations contribute less than 0.2% of the total operations in the and are easily implemented using low-cost DSP equipment, e.g., commensurate with current-generation Cortex M4 ARM processors. [0153] The total operations needed to detect and determine preliminary geo-observables using phase SCORE contribute less than 20% of the total operations in the processor, and are implementable using low-to-medium cost DSP equipment, e.g., Cortex M4 or A8 ARM processors, for most of the processing scenarios shown in this Table. [0154] All of the operations shown in these scenarios are implementable using currently available DSP equipment, e.g., Cortex A18 ARM processors. Moreover, although the complexity of conjugate SCORE can be high, it scales linearly with the number of emitters detected in using phase-SCORE, and is dominated by highly regular operations, such as the CCM computation, that can be implemented in FPGA or ASIC's.

[0155] It should be noted as well that many of these operations, or operations of similar complexity, are performed in the general civil GNSS flow diagram shown in FIG. 20. In this regard, the summary complexity through phase-SCORE detection operations is likely to be commensurate with processing requirements for any civil GNSS signals employing unmodulated pilot signals, or in resilient PNT use scenario's in which baseband symbol data can be supplied to the processor, e.g., to develop PNT analytics in order to aid a primary PNT unit.

DETAILED DESCRIPTION OF VARIOUS ASPECTS OF THE INVENTION

Baseline Single-Antenna Reception

[0156] The Blind Despreading Anti-Spoofing (BDAS) approach for the baseline single-feed reception scenario is shown in FIG. 1, in which a GNSS receiver (1) with position modeled by 3×1 time-varying vector p.sub.R(t) and yaw-pitch-roll orientation modeled by 3×3 time-varying rotation matrix Ψ.sub.R(l) is assumed to be operating within the field of view (FoV) of L.sub.T transmitters (3,5,7), e.g., satellite vehicles (SV's) (3), with positions modeled by 3×1 time-varying vectors {p.sub.T(t; l)}.sub.l=1.sup.L.sup.T. Ignoring relativistic effects and atmospheric propagation factors, the observed position of transmitter within the frame of reference of the receiver can then be modeled by 3×1 position vector p.sub.TR(t; l)=p.sub.T(t; l)−p.sub.R(t). Over short time intervals, p.sub.TR(t; l) can be further approximated by

[00006]$\begin{array}{cc}{p}_{\mathrm{TR}}\left(t;\right)\approx p_{\mathrm{TR}}\left(t_0;\right)+v_{\mathrm{TR}}\left(t_0;\right).Math.\left(t-t_0\right)+\frac{1}{2}.Math.a_{\mathrm{TR}}\left(t_0;\right).Math.{\left(t-t_0\right)}^2,& \mathrm{Eqn}.Math..Math.\left(4\right)\end{array}$

where v.sub.TR(t; l)=v.sub.T(t; l)−v.sub.R(t) and a.sub.TR(t; l)=a.sub.T(t; l)−a.sub.R(t) are the observed velocity and acceleration of transmitter l within the frame of reference of the receiver, and where |t−t.sub.0| is small.

[0157] Each GNSS transmitter generates a short civil GNSS signal-in-space (SiS), given by √{square root over (2P.sub.T(l))}Re{s.sub.T(t−τ.sub.T(l); l)e.sup.j(2πf.sup.T.sup.t+φ.sup.T.sup.(l))} for transmitter l in FIG. 1, where f.sub.T(l) is the transmit frequency, typically common to all SV's (f.sub.T(l)≡f.sub.T), and where P.sub.T(l), φ.sub.T(l) and τ.sub.T(l) are the transmit power, phase, and electronics delay for SiS l, respectively.

[0158] As shown in FIG. 2, the civil GNSS complex baseband signal s.sub.T(t; l) is a modulation-on-symbol direct-sequence spread spectrum (MOS-DSSS) signal, given by

[00007]$\begin{array}{cc}{s}_T\left(t;\right)=\underset{n_{\mathrm{sym}}}{.Math.}.Math..Math.d_T\left(n_{\mathrm{sym}};\right).Math.h_T\left(t-T_{\mathrm{sym}}.Math.n_{\mathrm{sym}};\right),& \mathrm{Eqn}.Math..Math.\left(5\right)\end{array}$

where d.sub.T(n.sub.sym; l) is a civil GNSS baseband symbol stream with symbol rate f.sub.sym=1/T.sub.sym, and h.sub.T(t; l) is a spread symbol given by

[00008]$\begin{array}{cc}{h}_T\left(t;\right)=\underset{m_{\mathrm{chp}}=0}{\overset{M_{\mathrm{chp}}}{.Math.}}.Math..Math.c_T\left(m_{\mathrm{chp}};\right).Math.g_T\left(t-T_{\mathrm{chp}}.Math.m_{\mathrm{chp}};\right),.Math.T_{\mathrm{chp}}=\frac{T_{\mathrm{sym}}}{M_{\mathrm{chp}}},& \mathrm{Eqn}.Math..Math.\left(6\right)\end{array}$

and where {c.sub.T(m.sub.chp; l)}.sub.m.sub.chp.sub.=0.sup.M.sup.chp.sup.−1 is a pseudorandom code with length M.sub.chp, and g.sub.T(t; l) is the chip symbol modulated by each chip. For the GPS L1 signal, the chip shaping is nominally rectangular and identical for each transmitted signal (g.sub.T(t; l)≡g.sub.T(t)); however, in practice, the shaping can vary between transmitters. The chip shaping can also vary between different GNSS formats.

[0159] For the L1 GPS C/A signal, the symbol period is 1 millisecond (ms) and the baseband symbol rate is 1 kilosample/second (ksps), the ranging code is a Gold code with pseudorandom noise (PRN) index k.sub.PRN(l) and length M.sub.chp=1,023, and the chip-shaping is a rectangular pulse with duration T.sub.chp=1/1,023 ms, modulated by digital and analog processing elements in the transmission chain. The L1 GPS symbol stream also possess additional exploitable structure, as it is constructed from a 50 bit/second (bps) BPSK (real) navigation signal b.sub.T(n.sub.NAV; l), repeated 20 times per MOS-DSSS symbol to generated the 1 ksps baseband MOS-DSSS symbol sequence d.sub.T(n.sub.sym; l), i.e., d.sub.T(n.sub.sym; l)=b.sub.T(└n.sub.sym/M.sub.NAV┘; l) where M.sub.NAV=20 and └•┘ denotes the integer truncation operation. Moreover, the navigation signal possess internal fields that can be exploited to aid the demodulation process, e.g., the TLM Word Preamble transmitted every 6 seconds, i.e., within every 300-bit Navigation subframe.

[0160] This model ignores other SiS's that may be transmitted simultaneously with civil GNSS signals, e.g., the GPS P(Y), M, and L1C SiS's. In practice, these signals are largely outside the band occupied by the civil GNSS signals, and are moreover received well below the noise floor. Therefore, their effect on the despreading methods described herein should be small. An important exception is military GNSS signals that may be generated by ground pseudolites, as those signal components may be above the noise floor within the civil GNSS passband in some use scenarios.

[0161] A list of civil GNSS signals that possess MOS-DSSS modulation format is provided in FIG. 3. As this Figure shows, every GNSS network deployed to date possesses at least one signal type that can be modeled by Eqn (5)-Eqn (6). Moreover, all of these signals have a 1 millisecond (ms) ranging code period, albeit with widely varying chip rates and SiS bandwidths, allowing their demodulation using a common receiver structure. Lastly, as the last column of FIG. 3 shows, the symbol streams of each of these signals possess properties that can be used to detect and differentiate them from other signals in the environment, including military GNSS signals, jammers, and civil GNSS signals with the same or similar properties.

[0162] Using either of the single-feed receiver front-ends shown in FIG. 4 and FIG. 5 and described in detail earlier, or equivalent reception structures, the single-feed GNSS receiver (1) receives a SiS over a radio frequency (RF) band encompassing at least part of a civil GNSS band; downconverts that signal by frequency shift f.sub.R to complex baseband, where it is represented as a complex signal x.sub.R(t−τ.sub.R) comprising an in-phase (I) rail as its real part and a quadrature (Q) rail as its imaginary part; or downconverts that signal by frequency shift f.sub.R−f.sub.IF to an intermediate frequency (IF) f.sub.IF, where it is represented as a real signal; and where a carrier phase offset φ.sub.R and delay τ.sub.R are induced in the signal due to impairments in the receiver's local oscillator (LO) and reception electronics. The downconverted analog signal x.sub.R(t−τ.sub.R) is then converted to digital format using an analog-to-digital converter (ADC) with sampling rate M.sub.ADCf.sub.sym, where f.sub.sym is the civil GNSS signal symbol rate shown in FIG. 2 and M.sub.ADC is a positive integer.

[0163] The receiver frequency shift f.sub.R need not be equal to the civil GNSS transmit frequency f.sub.T (or f.sub.T(l)), and the receiver bandwidth may encompass only a portion of the civil GNSS band of interest to the receiver. Moreover, the ADC symbol multiplier M.sub.ADC need not bear any relationship to the spreading factor M.sub.chp shown in FIG. 2. This can greatly improve cost and flexibility of the receiver hardware. For example, if the receiver is implemented using a software-defined radio (SDR) architecture as shown in FIG. 4B or FIG. 5B, the receiver can be flexibly tasked to any civil GNSS band, or to any portion of a civil GNSS band, and can downconvert and sample that band, or that portion of a band, using exactly the same hardware, regardless of the parameters of the civil GNSS signal being transmitted on that band. Moreover, the receiver can use low-cost, commercial off-the-shelf components with generic bandwidths, downconversion frequencies, and sampling rates, rather than custom built gear that is matched to the transmit frequency, bandwidth, and chip-rate of the civil GNSS signal.

[0164] Assuming the direct-conversion receiver structure shown in FIG. 4, x.sub.R(l) is given by

[00009]$\begin{array}{cc}{x}_R\left(t\right)=i_R\left(t\right)+\underset{=1}{\overset{L}{.Math.}}.Math..Math.s_R\left(t;\right),& \mathrm{Eqn}.Math..Math.\left(7\right)\end{array}$

where i.sub.R(t) is the noise and interference received in the receiver passband and s.sub.R(t; l) is the complex-baseband representation of the signal received from transmitter l. In absence of multipath, and ignoring atmospheric propagation effects, s.sub.R(t; l) is modeled by

s.sub.R(t;l)=g.sub.TR(t;l)e.sup.j2π(f.sup.T.sup.(l)−f.sup.R.sup.)te.sup.−j2πf.sup.T.sup.(l)τ.sup.TR.sup.(t;l)s.sub.T(t−τ.sub.TR(t;l)−τ.sub.T(l);l),   Eqn. (8)

g.sub.Tr(t;l)=√{square root over (P.sub.R(t;l))}e.sup.j(φ.sup.T.sup.(l)−φ.sup.R.sup.),   Eqn. (9)

where

[00010]${\tau }_{\mathrm{TR}}\left(t;\right)=\frac{1}{c}.Math.{.Math.p_{\mathrm{TR}}\left(t;\right).Math.}_2$

and P.sub.R(t; l) are the observed time-of-arrival (TOA) (defined separately from the aggregate delay τ.sub.T(l) induced in the transmitter electronics) and received incident power (RIP) of GNSS signal l at the receiver, respectively, and where ∥•∥.sub.2 denotes the Euclidean L-2 norm. The time-varying TOA also induces a Doppler shift in s.sub.R(t; l), given by frequency of arrival (FOA)

α.sub.TR(t;l)=−f.sub.T(l)τ.sub.TR.sup.(1)(t;l),   Eqn. (10)

defined separately from the frequency offset f.sub.T(l)−f.sub.R induced by difference between the transmit and target receive frequency, where τ.sub.TR.sup.(1)(t; l) is the first derivative of τ.sub.TR(t; l),

[00011]$\begin{array}{cc}{\tau }_{\mathrm{TR}}^{\left(1\right)}\left(t;\right)=\frac{1}{c}.Math.\frac{d}{\mathrm{dt}}.Math.{.Math.p_{\mathrm{TR}}\left(t;\right).Math.}_2=\frac{1}{c}.Math.u_{\mathrm{TR}}^T\left(t;\right).Math.v_{\mathrm{TR}}\left(t;\right),& \mathrm{Eq}.Math..Math.\left(11\right)\end{array}$

where u.sub.TR(t; l)=p.sub.TR(t; l)/∥p.sub.TR(t; l∥.sub.2 is the line-of-bearing (LOB) direction vector from the receiver to transmitter l, which can be expressed in azimuth relative to East and elevation relative to the plane of the Earth using well-known directional transformations. The local direction-of-arrival (DOA) of the signal received from transmitter l is then represented by local DOA direction vector u.sub.R(t; l)=Ψ.sub.R(t)u.sub.TR(t; l), which can be converted to azimuth relative to the local receiver heading and elevation relative to the plane of receiver motion using well-known directional transformations. Lastly, if the transmitter and receiver phases (assumed constant over the reception interval) can be calibrated and removed from the end-to-end phase measurement, then the calibrated carrier phase (CCP) φ.sub.TR(t; l)=−2πf.sub.T(l)τ.sub.TR(t; l) can be used to provide a precise estimate of range from the transmitter to the receiver.

[0165] The TOA, RIP, FOA, LOB, DOA, and CCP are summarized in FIG. 6, where they are referred to as geo-observable, as they are observable parameters of the GNSS received signals that can be used to geo-locate the GNSS receiver, given knowledge of the GNSS transmitter locations over time, which are transmitted within the navigation stream of those signals. They are also referred to here as pointing, navigation, and timing (PNT) analytics, as they provide metadata of the transmitted signals that can be used for purposes beyond geo-location of the receiver, e.g., to assess quality and reliability of the received GNSS signals, in order to aid other PNT systems on board the receiver platform.

[0166] In this regard, the following reception scenarios are particularly pertinent to the invention: [0167] SDR reception scenarios, in which the invention is used to rapidly scan all of the GNSS bands, and assess availability of GNSS resources for dedicated GNSS receivers already on board the receiver platform. A feature of the invention that is particularly useful in this scenario is its ability to provide a fast and robust time-to-first-fix (TTFF) from a cold start, through its ability to detect all of the GNSS transmitters in the receiver's field of view (FoV), and to determine usable geo-observables of signals received from those transmitters, without a search over timing or ranging codes, using a flexible SDR receiver with a standardized reception bandwidth. [0168] Beacon reception scenarios, in which the invention is used to separate signals from terrestrial beacons. This can include mixed GNSS SV, beacon reception scenarios, in which the terrestrial beacons are operating in the same band as the GNSS SV's, as well as scenarios in which the beacons are operating in separate bands. A feature of the invention that is particularly useful in this scenario is its ability to detect, separate, and determine usable geo-observables of GNSS signals with extreme wide disparity in RIP, and with receive signal-to-interference ratios that are strong enough to substantively degrade performance of conventional correlative signal detection methods. [0169] Targeted spoofing scenarios, in which the invention is used to separate signals from GNSS SV's and spoofers that are attempting to deceive the receiver, by generating signals that jam the legitimate GNSS SV's and emulate those SV's (aligned spoofers) or SV's in other geographical locations. A feature of the invention that is particularly useful in this scenario is its ability to detect and separate GNSS signals with close (covert spoofing) or wide (non-covert spoofing) disparity in RIP, and to estimate geo-observables that can be used to both detect spoofing, and (with appropriate ancillary information) identify which signals are emanating from spoofers or legitimate GNSS transmitters.

[0170] Detection and geo-observable parameter estimation methods developed in the context of these applications are referred to here as resilient PNT analytics.

[0171] As shown in FIG. 10 for the dynamic mixed GNSS SV/beacon reception environment shown in FIG. 9, although the geo-observables are in general time-varying, they adhere closely to a zero-order model in RIP, and to a first-order model in TOA and FOA, over time intervals on the order of 1-10 seconds, such that

P.sub.R(t;l)=P.sub.R(t.sub.0;l),   Eqn (12)

τ.sub.TR(t;l)≈τ.sub.TR(t.sub.0;l)+τ.sub.TR.sup.(1)(t.sub.0;l)(t−t.sub.0),   Eqn. (13)

α.sub.TR(t;l)≈α.sub.TR(t.sub.0;l)+α.sub.TR.sup.(1)(t.sub.0;l)(t−t.sub.0),   Eqn. (14)

for t.sub.0≦t≦t.sub.0+10 s. Taking all of the phase and delay shifts induced in the transmitter and receiver hardware together yields received link l GNSS signal model

[00012]$\begin{array}{cc}{s}_R\left(t_0+t;l\right)\approx e^{j.Math..Math.2.Math..Math.\pi \left({\alpha }_{\mathrm{TR}}\left(l\right).Math.t+\frac{1}{2}.Math.{\alpha }_{\mathrm{TR}}^{\left(1\right)}\left(l\right).Math.t^2\right)}.Math.g_{\mathrm{TR}}\left(l\right).Math.s_T\left(\left(1-{\tau }_{\mathrm{TR}}^{\left(1\right)}\left(l\right)\right).Math.t-{\tau }_{\mathrm{TR}}\left(l\right);l\right),& \mathrm{Eqn}.Math..Math.\left(15\right)\end{array}$

at the input to the dual-ADC's shown in FIG. 4, where

g.sub.TR(l)√{square root over (P.sub.R(t.sub.0;l))}exp(jφ.sub.TR(l)), end-to-end complex channel gain   Eqn (16)

φ.sub.TR(l)φ.sub.T(l)−φ.sub.R−2πf.sub.T(l)τ.sub.TR(t.sub.0;l), end-to-end channel phase   Eqn (17)

τ.sub.TR(l)τ.sub.T(l)+τ.sub.R+τ.sub.TR(t.sub.0;l), end-to-end observed TOA   Eqn (18)

α.sub.TR(l)f.sub.T(l)+f.sub.R+α.sub.TR(t.sub.0;l), end-to-end observed FOA   Eqn (19)

τ.sub.TR.sup.(1)(l)τ.sub.TR.sup.(1)(t.sub.0;l), observed differential TOA (DTOA)   Eqn (20)

α.sub.TR.sup.(1)(l)α.sub.TR.sup.(1)(t.sub.0;l), observed differential FOA (DFOA).   Eqn (21)

Similarly, GNSS signal l can be approximated by

{tilde over (s)}.sub.R(t.sub.0+t;l)≈√{square root over (2)}Re{e.sup.j2πf.sup.IF.sup.ts.sub.R(t.sub.0+t;l)},   Eqn (22)

at the input to the ADC shown in FIG. 5, where s.sub.R(t.sub.0+t; l) is given by Eqn (15) and φ.sub.R=φ.sub.R.sup.(1)+φ.sub.R.sup.(2).

[0172] As shown in FIG. 11, and FIG. 12 the LOB and local DOA also adhere closely to a zero-order and first-order model over similar time intervals. Moreover, the variation in local DOA is nearly identical for all of the received GNSS signals, as it primarily due to motion/orientation of the receiver, which equally affects all of the signal LOB's.

[0173] As shown in FIG. 13, the ADC output signal is passed to a symbol-rate synchronous vector channelizer, integral to the resent invention, which transforms that signal into a M.sub.DoF×1 vector signal x.sub.DoF(n.sub.sym) with rate f.sub.sym equal to the rate of the baseband symbol sequence of the MOS-DSSS civil GNSS signal shown in FIG. 2, e.g., 1 ksps for the civil GNSS signals shown in FIG. 3, and where M.sub.DoF is the degrees of freedom (DoF's) of the channelizer output signal. The symbol-rate synchronous channelizer shown in FIG. 13 is implemented using a 1:M.sub.ADC serial-to-parallel (S/P) operation, which forms M.sub.ADC×1 S/P output signal x.sub.sym(n.sub.sym) given by Eqn (1), where T.sub.sym=1/f.sub.sym is the baseband symbol period, followed by multiplication by M.sub.DoF×M.sub.ADC channelization matrix T.sub.DoF, forming M.sub.DoF×1 channelizer output signal x.sub.DoF(n.sub.sym)=T.sub.DoF×sym(n.sub.sym). Both signals have data rate f.sub.sym.

[0174] FIG. 14 describes a particularly useful channelizer based on fast Fourier transform (FFT) operations. The channelizer first applies FFT with M.sub.ADC input samples, M.sub.FFT≧M.sub.ADC output bins, and a (nominally) rectangular) window {w.sub.FFT(m.sub.ADC)}.sub.m.sub.ADC.sub.=0.sup.M.sup.ADC.sup.−1 to each S/P output vector. The channelizer then selects M.sub.DoF output bins {k.sub.bin(m.sub.DoF)}.sub.m.sub.DoF.sub.=0.sup.M.sup.DoF.sup.−1 for use in subsequent adaptive processing algorithms. The channelizer thereby implements channelizer matrix

[00013]$T_{\mathrm{DoF}}=\left[\exp \left(-j.Math..Math.2.Math.\pi .Math..Math.k_{\mathrm{bin}}\left(m_{\mathrm{DoF}}\right).Math.\frac{m_{\mathrm{ADC}}}{M_{\mathrm{FFT}}}\right).Math.w_{\mathrm{FFT}}\left(m_{\mathrm{ADC}}\right)\right].$

The output bins can be preselected, e.g., to avoid known receiver impairments such as LO leakage, or to reduce complexity of the channelization operation, e.g., using sparse FFT methods; or adaptively selected, e.g., to avoid dynamic narrowband co-channel interferers (NBCCI).

[0175] In one embodiment, the approach is extended polyphase filtering operation x.sub.DoF(n.sub.sym)=

[00014]$T_{\mathrm{DoF}}.Math.\underset{m_{\mathrm{sym}}=0}{\overset{M_{\mathrm{sym}}}{.Math.}}.Math..Math.\left(w_{\mathrm{sym}}\left(m_{\mathrm{sym}}\right)\odot x_{\mathrm{sym}}\left(n_{\mathrm{sym}}-m_{\mathrm{sym}}\right)\right),$

where {w.sub.sym(w.sub.sym(m.sub.sym)}.sub.m.sub.sym.sub.=0.sup.M.sup.sym.sup.−1 is an order-M.sub.sym polyphase filter, “⊙” denotes the element-wise multiplication operation, and M.sub.DoF×M.sub.ADC matrix

[00015]$T_{\mathrm{DoF}}=\left[\exp \left(-j.Math..Math.2.Math.\pi .Math..Math.k_{\mathrm{bin}}\left(m_{\mathrm{DoF}}\right).Math.\frac{m_{\mathrm{ADC}}}{M_{\mathrm{FFT}}}\right)\right]$

implements an unpadded, unwindowed FFT with M.sub.ADC input samples and M.sub.DoF output bins. This channelizer is especially useful in environments subject to very strong NBCCI. However, it affects dispersion added to the channelizer output signal.

[0176] Given the reception scenario described in FIG. 6, then the single-sample channelizer output signal shown in FIG. 13 can be modeled by

[00016]$\begin{array}{cc}{x}_{\mathrm{DoF}}\left(n_{\mathrm{sym}}\right)=i_{\mathrm{DoF}}\left(n_{\mathrm{sym}}\right)+\underset{=1}{\overset{L_T}{.Math.}}.Math..Math.s_{\mathrm{DoF}}\left(n_{\mathrm{sym}};\right),& \mathrm{Eqn}.Math..Math.\left(23\right)\end{array}$

where i.sub.DoF(n.sub.sym)=T.sub.DoFi.sub.sym(n.sub.sym) and s.sub.DoF(n.sub.sym; l)=T.sub.DoFs.sub.sym(n.sub.sym; l) are the M.sub.DoF×1 interference and GNSS signal l channelizer output signals, respectively, and where

[00017]$i_{\mathrm{sym}}\left(n_{\mathrm{sym}}\right)=\left[i_R\left(T_{\mathrm{sym}}\left(n_{\mathrm{sym}}+\frac{m_{\mathrm{ADC}}}{M_{\mathrm{ADC}}}\right)\right)\right].Math..Math.\mathrm{and}$$s_{\mathrm{sym}}\left(n_{\mathrm{sym}};\right)=\left[s_R\left(T_{\mathrm{sym}}\left(n_{\mathrm{sym}}+\frac{m_{\mathrm{ADC}}}{M_{\mathrm{ADC}}}\right);\right)\right]$

are the M.sub.ADC×1 interference and GNSS signal l S/P output vectors, respectively. Assuming the RIP, TOA, and FOA time variation models given in Eqn (12)-Eqn (14), decomposing τ.sub.TR(l) and α.sub.TR(l) into symbol-normalized FOA and TOA components

[00018]$\begin{array}{cc}{\tau }_{\mathrm{TR}}\left(\right)=\left({\tilde{\tau }}_{\mathrm{TR}}\left(\right)+n_{\mathrm{TR}}\left(\right)\right).Math.T_{\mathrm{sym}}.Math.\{\begin{array}{cc}{\tilde{\tau }}_{\mathrm{TR}}\left(\right)\in \left[0.Math..Math.1\right),& \mathrm{fractional}.Math..Math.\mathrm{fine}.Math..Math.\mathrm{initial}.Math..Math.\mathrm{TOA}\\ n_{\mathrm{TR}}\left(\right)\in \mathbb{Z},& \mathrm{initial}.Math..Math.\mathrm{TOA}.Math..Math.\mathrm{symbol}.Math..Math.\mathrm{offset},\end{array}& \mathrm{Eqn}.Math..Math.\left(24\right)\\ {\alpha }_{\mathrm{TR}}\left(\right)=\left({\tilde{\alpha }}_{\mathrm{TR}}\left(\right)+k_{\mathrm{TR}}\left(\right)\right).Math.f_{\mathrm{sym}}.Math.\{\begin{array}{cc}{\tilde{\alpha }}_{\mathrm{TR}}\left(\right)\in \left[-\frac{1}{2}.Math..Math.\frac{1}{2}\right),& \mathrm{fractional}.Math..Math.\mathrm{fine}.Math..Math.\mathrm{FOA}\\ k_{\mathrm{TR}}\left(\right)\in \mathbb{Z},& \mathrm{FOA}.Math..Math.\mathrm{Nyquist}.Math..Math.\mathrm{zone},\end{array}& \mathrm{Eqn}.Math..Math.\left(25\right)\end{array}$

and defining dimensionless symbol-normalized DFOA {tilde over (α)}.sub.TR.sup.(l)(l) α.sub.TR.sup.(1)(l)T.sub.sym.sup.2, then over short (<10 s) reception intervals s.sub.DoF(n.sub.sym; l) can be modeled by

[00019]$\begin{array}{cc}{s}_{\mathrm{DoF}}\left(n_{\mathrm{sym}};\right)=A_{\mathrm{DoF}}\left(n_{\mathrm{sym}};\right).Math.d_T\left(n_{\mathrm{sym}}-n_{\mathrm{TR}}\left(\right);\right).Math.e^{j.Math..Math.2.Math.\pi }\left({\tilde{\alpha }}_{\mathrm{TR}}\left(\right).Math.n_{\mathrm{sym}}+\frac{1}{2}.Math.{\tilde{\alpha }}_{\mathrm{TR}}^{\left(1\right)}\left(\right).Math.n_{\mathrm{sym}}^2\right),& \mathrm{Eqn}.Math..Math.\left(26\right)\end{array}$

where d.sub.T(n.sub.sym; l) is a Q.sub.sym(l)×1 transmitted symbol delay-line vector given by

[00020]$\begin{array}{cc}{d}_T\left(n_{\mathrm{sym}};\right)=\left(\begin{array}{c}{d}_T\left(n_{\mathrm{sym}}-Q_{\min }\left(\right);\right)\\ \mathrm{.Math.}\\ d_T\left(n_{\mathrm{sym}}-Q_{\max }\left(\right);\right)\end{array}\right),& \mathrm{Eqn}.Math..Math.\left(27\right)\end{array}$

such that {d.sub.T(n.sub.sym−q.sub.sym−n.sub.TR(l); l)}.sub.q.sub.sym.sub.=Q.sub.min.sub.(l).sup.Q.sup.max.sup.(l) are all the transmitted data symbols observable within S.sub.DoF(n.sub.sym; l), and A.sub.DoF(n.sub.sym; l) is a slowly time-varying, dispersive M.sub.DoF×Q.sub.sym(l) multiple-input, multiple-output (MIMO) link signature matrix between the transmitted GNSS l symbol stream and the channelizer output vector, given by

[00021]$\begin{array}{cc}{A}_{\mathrm{DoF}}\left(n_{\mathrm{sym}};\right)=T_{\mathrm{DoF}}.Math.A_{\mathrm{sym}}\left(n_{\mathrm{sym}};\right),& \mathrm{Eqn}.Math..Math.\left(28\right)\\ A_{\mathrm{sym}}\left(n_{\mathrm{sym}};\right)=g_{\mathrm{TR}}\left(\right).Math.{\Delta }_R\left(n_{\mathrm{sym}};\right).Math.H_R\left(n_{\mathrm{sym}};\right),& \mathrm{Eqn}.Math..Math.\left(29\right)\\ {\Delta }_R\left(n_{\mathrm{sym}};\right)=\mathrm{diag}.Math.\left\{\exp (.Math.j.Math..Math.2.Math.\pi \left({\tilde{\alpha }}_{\mathrm{TR}}\left(\right)+k_{\mathrm{TR}}\left(\right)+{\tilde{\alpha }}_{\mathrm{TR}}^{\left(1\right)}\left(\right).Math.\left(n_{\mathrm{sym}}+\frac{1}{2}.Math.\frac{m_{\mathrm{ADC}}}{M_{\mathrm{ADC}}}\right)\right).Math.\frac{m_{\mathrm{ADC}}}{M_{\mathrm{ADC}}})\right\},& \mathrm{Eqn}.Math..Math.\left(30\right)\\ H_R\left(n_{\mathrm{sym}};\right)=\hspace{1em}\left[h_T\left(t_0+T_{\mathrm{sym}}\left(q_{\mathrm{sym}}+\frac{m_{\mathrm{ADC}}}{M_{\mathrm{ADC}}}-{\tilde{\tau }}_{\mathrm{TR}}\left(\right)-{\tau }_{\mathrm{TR}}^{\left(1\right)}\left(\right).Math.\left(n_{\mathrm{sym}}+\frac{m_{\mathrm{ADC}}}{M_{\mathrm{ADC}}}\right)\right);\right)\right],& \mathrm{Eqn}.Math..Math.\left(31\right)\end{array}$

and where Q.sub.sym(l)=Q.sub.max(l)−Q.sub.min(l)+1. Note that Eqn (30) is a function of the symbol-normalized normalized DFOA {tilde over (α)}.sub.TR.sup.(1)(l) and fine fractional FOA {tilde over (α)}.sub.TR(l), while Eqn (31) is a function of the DTOA τ.sub.TR.sup.(1)(l) and symbol-normalized TOA {tilde over (τ)}.sub.TR(l).

[0177] Over short observation intervals, and assuming the single-symbol channelization approach shown in FIG. 13, Q.sub.min(l)=0, and Q.sub.max(l)=1 or 2, resulting in Q.sub.sym(l)=2 or 3, with the larger value only occurring if {tilde over (τ)}.sub.TR(l)≈1. Over sufficiently long observation intervals, Q.sub.max(l).fwdarw.+2 if ε.sub.TR(l)>0, or Q.sub.min(l).fwdarw.−1 if ε.sub.TR(l)<0. Thus all reasonable TOA drift can be accommodated over long observation intervals by setting Q.sub.min(l)≡−1 and Q.sub.max(l)≡+2, yielding Q.sub.sym(l)≡4, for all of the GNSS signals in the receiver's field of view. However, even in this case, the lagging or leading column of A.sub.DoF(n.sub.sym; l) is substantively weaker than the base columns, such that A.sub.DoF(n.sub.sym; l) has “substantive” rank 2. If the channelizer spans multiple symbols, e.g., if it is implemented using an order-M.sub.sym polyphase filter, then Q.sub.max(l) will be increased by the number of additional symbols spanned in the channelizer, e.g., Q.sub.max(l)←Q.sub.max(l)+M.sub.sym and Q.sub.sym(l)←Q.sub.sym(l)+M.sub.sym.

[0178] Defining observed GNSS signal l received symbol vector

[00022]$\begin{array}{cc}{d}_R\left(n_{\mathrm{sym}};l\right).Math.\stackrel{\Delta }{=}.Math.d_T\left(n_{\mathrm{sym}}-n_{\mathrm{TR}}\left(l\right);l\right).Math.\exp \left(j.Math..Math.2.Math..Math.\pi \left({\tilde{\alpha }}_{\mathrm{TR}}\left(l\right).Math.n_{\mathrm{sym}}+\frac{1}{2}.Math.{\tilde{\alpha }}_{\mathrm{TR}}^{\left(1\right)}\left(l\right).Math.n_{\mathrm{sym}}^2\right)\right),& \mathrm{Eqn}.Math..Math.\left(32\right)\end{array}$

and substituting Eqn (26) into Eqn (23), yields channelizer output signal model

[00023]$\begin{array}{cc}\begin{array}{c}{x}_{\mathrm{DoF}}\left(n_{\mathrm{sym}}\right)=.Math.i_{\mathrm{DoF}}\left(n_{\mathrm{sym}}\right)+\underset{=1}{\overset{L_T}{.Math.}}.Math.A_{\mathrm{DoF}}\left(n_{\mathrm{sym}};\right).Math.d_R\left(n_{\mathrm{sym}};\right)\\ =.Math.i_{\mathrm{DoF}}\left(n_{\mathrm{sym}}\right)+A_{\mathrm{DoF}}\left(n_{\mathrm{sym}}\right).Math.d_R\left(n_{\mathrm{sym}}\right),\end{array}& \mathrm{Eqn}.Math..Math.\left(33\right)\end{array}$

where d.sub.R(n.sub.sym) is the Q.sub.net×1 dispersive network received symbol vector and A.sub.DoF(n.sub.sym) is the slowly time-varying, dispersive M.sub.chn×Q.sub.net MIMO network signature matrix between all of the GNSS transmitted and the received symbols, given by

[00024]$\begin{array}{cc}{A}_{\mathrm{DoF}}\left(n_{\mathrm{sym}}\right)=\left[A_{\mathrm{DOF}}\left(n_{\mathrm{sym}};1\right).Math..Math.\mathrm{.Math.}.Math..Math.A_{\mathrm{DOF}}\left(n_{\mathrm{sym}};L_T\right)\right],& \mathrm{Eqn}.Math..Math.\left(34\right)\\ d_R\left(n_{\mathrm{sym}}\right)=\left(\begin{array}{c}{d}_R\left(n_{\mathrm{sym}};1\right)\\ \mathrm{.Math.}\\ d_R\left(n_{\mathrm{sym}};L_T\right)\end{array}\right),& \mathrm{Eqn}.Math..Math.\left(35\right)\end{array}$

respectively, and where

[00025]$Q_{\mathrm{net}}.Math.\stackrel{\Delta }{=}.Math.\underset{=1}{\overset{L_T}{.Math.}}.Math.Q_{\mathrm{sym}}\left(\right)$

is the total network loading induced by the GNSS signals.

[0179] If M.sub.DoF≧Q.sub.net and the GNSS signals are received with high despread SINR, the entire network symbol stream can be extracted from the channel, i.e., the interfering GNSS symbol streams can be excised from each GNSS symbol, using purely linear combining operations. This can include methods well known to the signal processing community, e.g., blind adaptive baseband extraction methods described in the prior art, and linear minimum-mean-square-error (LMMSE) methods described in the prior art. For signals transmitted from MEO GNSS SV's, and in the absence of pseudolites or spoofers operating in the same band as those SV's, no more than 12 SV's are likely to be within the field of view of a GNSS receiver at any one time (L.sub.T≦12). Over short observation intervals, Q.sub.sym(l)≈2, and the maximum substantive rank of the MIMO network signature matrix is therefore likely to be 24. This number can grow to 48 in the presence of 12 spoofers “assigned” to each legitimate GNSS signal. In both cases, Q.sub.net is much less than the number of channels available for any GNSS signals listed in FIG. 3.

[0180] This channel response is also exactly analogous to channel responses induced in massive MIMO networks currently under investigation for next-generation (5G) cellular communication systems. However, it achieves this response using only a single-feed receiver front-end, thereby bypassing the most challenging aspect of massive MIMO transceiver technology. And it provides an output signal with an effective data rate of 1 ksps for all of the signals listed in FIG. 3, as opposed to 10-20 Msps for massive MIMO networks under consideration for 5G cellular communication applications. This receiver front-end also provides considerable flexibility in the number of degrees of freedom M.sub.DoF used at the receiver, e.g., by modifying the ADC sampling rate, or the number of FFT bins used in the channelizer, or any combination thereof (albeit, at loss in processing gain if the channelizer does not fully cover the GNSS signal passband).

[0181] Lastly, the digital signal processing (DSP) operations needed to exploit this channel response are expected to be very similar to operations needed for 5G data reception, albeit at a 3-4 order-of-magnitude lower switching rate. Given the massive investment expected in 5G communications over the next decade, and the ongoing exponential improvements in cost and performance of DSP processing and memory, e.g., Moore's and Kryder's Laws, the ability to fully exploit this channel response will become increasingly easier over time.

Modeling Link Signature Blur

[0182] Over long observation intervals, the time-varying components of Eqn (30) and Eqn (31) results in “FOA blur” and “TOA blur” of the link signature matrix, exactly analogous to “signature blur” caused by change in observed direction-of-arrival (DOA) of moving platforms in conventional MIMO networks. Over moderate observation intervals, e.g., 0≦t<T.sub.symN.sub.sym where {tilde over (α)}.sub.TR.sup.(1)(l)N.sub.sym<1 in Eqn (30) and τ.sub.TR.sup.(1)(l)T.sub.symN.sub.sym<1 in Eqn (31), this can be quantified by expanding Δ.sub.R(n.sub.sym; l) and H.sub.R(n.sub.sym; l) in Taylor-series expansions

[00026]$\begin{array}{cc}\begin{array}{c}{\Delta }_R\left(n_{\mathrm{sym}};\right)=.Math.\begin{array}{c}.Math.\left(\underset{q=0}{\overset{\infty }{.Math.}}.Math.\frac{1}{q!}.Math.\mathrm{diag}.Math.\left\{{\left(j.Math..Math.2.Math.\pi .Math.{\tilde{\alpha }}_{\mathrm{TR}}^{\left(1\right)}\left(\right).Math.\left(n_{\mathrm{sym}}+\frac{1}{2}.Math.\frac{m_{\mathrm{ADC}}}{M_{\mathrm{ADC}}}\right).Math.\frac{m_{\mathrm{ADC}}}{M_{\mathrm{ADC}}}\right)}^q\right\}\right)\\ {\Delta }_R\left(\right)\end{array}\\ \approx .Math.\left(\underset{q=0}{\overset{\infty }{.Math.}}.Math.\frac{1}{q!}.Math.{\left(j.Math..Math.2.Math.\pi .Math.{\tilde{\alpha }}_{\mathrm{TR}}^{\left(1\right)}\left(\right).Math.n_{\mathrm{sym}}\right)}^q.Math.\mathrm{diag}.Math.\left\{{\left(\frac{m_{\mathrm{ADC}}}{M_{\mathrm{ADC}}}\right)}^q\right\}\right).Math.{\Delta }_R\left(\right),\end{array}& \mathrm{Eqn}.Math..Math.\left(36\right)\\ {\Delta }_R\left(\right)=\mathrm{diag}.Math.\left\{\exp \left(j.Math..Math.2.Math.\pi \left({\tilde{\alpha }}_{\mathrm{TR}}\left(\right)+k_{\mathrm{TR}}\left(\right)\right).Math.\frac{m_{\mathrm{ADC}}}{M_{\mathrm{ADC}}}\right)\right\},& \mathrm{Eqn}.Math..Math.\left(37\right)\\ \begin{array}{c}{H}_R\left(n_{\mathrm{sym}};\right)=.Math.\underset{q=0}{\overset{\infty }{.Math.}}.Math.\frac{1}{q!}.Math.\mathrm{diag}.Math.\left\{{\left(-{\tau }_{\mathrm{TR}}^{\left(1\right)}\left(\right).Math.T_{\mathrm{sym}}\left(n_{\mathrm{sym}}+\frac{m_{\mathrm{ADC}}}{M_{\mathrm{ADC}}}\right)\right)}^q\right\}.Math.H_R^{\left(q\right)}\left(\right)\\ \approx .Math.\underset{q=0}{\overset{\infty }{.Math.}}.Math.\frac{1}{q!}.Math.{\left(-{\tau }_{\mathrm{TR}}^{\left(1\right)}\left(\right).Math.T_{\mathrm{sym}}.Math.n_{\mathrm{sym}}\right)}^q.Math.H_R^{\left(q\right)}\left(\right),\end{array}& \mathrm{Eqn}.Math..Math.\left(38\right)\\ H_R^{\left(q\right)}\left(\right)=\left[h_T^{\left(q\right)}\left(t_0+T_{\mathrm{sym}}\left(q_{\mathrm{sym}}+\frac{m_{\mathrm{ADC}}}{M_{\mathrm{ADC}}}-{\tilde{\tau }}_{\mathrm{TR}}\left(\right)\right);\right)\right],& \mathrm{Eqn}.Math..Math.\left(39\right)\end{array}$

where h.sub.T.sup.(q)(t; l)=d.sup.qh.sub.T(t; l)/dt.sup.q. For the range of DTOA and DFOA values measured over the reception scenario shown in FIG. 9 and summarized in FIG. 10, |τ.sub.TR.sup.(1)(l)T.sub.sym| and |{tilde over (α)}.sub.TR.sup.(1)(l)| are less than 2.3×10.sup.−6 and 3.1×10.sup.−6, respectively, or 0.023 and 0.031 over 10,000 symbols (10 seconds), respectively, and Eqn (36)-Eqn (39) should therefore hold closely for small expansion orders. In this case, A.sub.DoF(n.sub.sym; l) is well represented by a low-order regression model over small observation intervals, e.g.,

[00027]$\begin{array}{cc}{A}_{\mathrm{DoF}}\left(n_{\mathrm{sym}};\right)=\underset{q_{\mathrm{blur}}=0}{\overset{Q_{\mathrm{blur}}\left(\right)}{.Math.}}.Math.A_{\mathrm{DoF}}^{\left(q_{\mathrm{blur}}\right)}\left(\right).Math.g_{\mathrm{blur}}^{\left(q_{\mathrm{blur}}\right)}\left(n_{\mathrm{sym}}\right),& \mathrm{Eqn}.Math..Math.\left(40\right)\end{array}$

where {g.sub.blur.sup.(q.sup.blur.sup.)(n.sub.sym)}.sub.q.sub.blur.sub.=0.sup.Q.sup.blur.sup.(l) is a deterministic sequence with order Q.sub.blur(l). Using Eqn (40), s.sub.DoF(n.sub.sym; l) is expressed as

[00028]$\begin{array}{cc}\begin{array}{c}{s}_{\mathrm{DoF}}\left(n_{\mathrm{sym}};\right)\approx .Math.\left(\underset{q_{\mathrm{blur}}=0}{\overset{Q_{\mathrm{blur}}\left(\right)}{.Math.}}.Math.A_{\mathrm{DoF}}^{\left(q_{\mathrm{blur}}\right)}\left(\right).Math.g_{\mathrm{blur}}^{\left(q_{\mathrm{blur}}\right)}\left(n_{\mathrm{sym}}\right)\right).Math.d_R\left(n_{\mathrm{sym}};\right)\\ =.Math.\begin{array}{c}\left[A_{\mathrm{DoF}}^{\left(0\right)}\left(\right).Math..Math.\mathrm{.Math.}.Math..Math.A_{\mathrm{DoF}}^{\left(Q_{\mathrm{blur}}\left(\right)\right)}\left(\right)\right]\\ \left(\begin{array}{c}{g}_{\mathrm{blur}}^{\left(0\right)}\left(n_{\mathrm{sym}}\right).Math.d_R\left(n_{\mathrm{sym}};\right)\\ \mathrm{.Math.}\\ g_{\mathrm{blur}}^{\left(Q_{\mathrm{blur}}\left(\right)\right)}\left(n_{\mathrm{sym}}\right).Math.d_R\left(n_{\mathrm{sym}};\right)\end{array}\right)\end{array}\\ =.Math.\begin{array}{c}\left[A_{\mathrm{DoF}}^{\left(0\right)}\left(\right).Math..Math.\mathrm{.Math.}.Math..Math.A_{\mathrm{DoF}}^{\left(Q_{\mathrm{blur}}\left(\right)\right)}\left(\right)\right]\\ \left(g_{\mathrm{blur}}.Math.\left(n_{\mathrm{sym}};\right).Math.d_R\left(n_{\mathrm{sym}};\right)\right),\end{array}\end{array}& \mathrm{Eqn}.Math..Math.\left(41\right)\end{array}$

where g.sub.blur(n.sub.sym; l)=[g.sub.blur.sup.(q.sup.blur)(n.sub.sym)].sub.q.sub.blur.sub.=0.sup.Q.sup.blur.sup.(l) is a (1+Q.sub.blur(l))×1 deterministic basis function and “” denotes the Kronecker product operation. A more useful form of Eqn (41) can be developed by reversing the order of the Kronecker product order, yielding link model

s.sub.DoF(n.sub.sym;l)≈A′.sub.DoF(l)d′.sub.R(n.sub.sym;l),   Eqn (42)

d′.sub.R(n.sub.sym;l)=d.sub.R(n.sub.sym;l)g.sub.blur(n.sub.sym;l),   Eqn (43)

A′.sub.DoF(l)=[A′.sub.DoF(Q.sub.min(l);l) . . . A′.sub.DoF(Q.sub.max(l);l)],   Eqn (44)

where A′.sub.DoF(q; l)=[a.sub.DoF.sup.(0)(:,q;l) . . . a.sub.DoF.sup.(Q.sup.blur.sup.(l))(:,q;l)] is the M.sub.DoF×(1+Q.sub.blur(l)) blurred signature matrix multiplying g.sub.blur(n.sub.sym; l)d.sub.R(n.sub.sym−q; l). Substituting Eqn. (42) into Eqn (33) yields network model

[00029]$\begin{array}{cc}\begin{array}{c}{x}_{\mathrm{DoF}}\left(n_{\mathrm{sym}}\right)=.Math.i_{\mathrm{DoF}}\left(n_{\mathrm{sym}}\right)+\underset{=1}{\overset{L_T}{.Math.}}.Math.A_{\mathrm{DoF}}^{\prime }\left(\right).Math.d_R^{\prime }\left(n_{\mathrm{sym}};\right)\\ =.Math.i_{\mathrm{DoF}}\left(n_{\mathrm{sym}}\right)+A_{\mathrm{DoF}}^{\prime }.Math.d_R^{\prime }\left(n_{\mathrm{sym}}\right).\end{array}& \mathrm{Eqn}.Math..Math.\left(45\right)\end{array}$

[0183] The rank of the MIMO link and network matrices grow to Q′.sub.sym(l)=Q.sub.sym(l)(1+Q.sub.blur(l)) and

[00030]$Q_{\mathrm{net}}^{\prime }=\underset{=1}{\overset{L_T}{.Math.}}.Math.Q_{\mathrm{sym}}^{\prime }\left(\right),$

respectively. Regression analysis of A.sub.DoF(n.sub.sym; l) for GPS C/A ranging codes and a channelizer output bandwidth of 250 kHz (M.sub.DoF˜250) shows that A.sub.DoF(n.sub.sym; l) is well modeled by Eqn (40) over observation intervals on the order of 2 seconds using a 1.sup.st-order polynomial sequence, i.e., Q.sub.blur(l)=1, if the signal is transmitted from a MEO SV, and using a 0.sup.th-order polynomial sequence, i.e., Q.sub.blur(l)=0, if the receiver is fixed and the symbol stream is transmitted from a pseudolite or ground beacon. Thus over long observation intervals the rank of the MIMO network matrix can double, e.g., to 48 in non-spoofing scenarios, or 96 in spoofing scenarios, which is still much lower than channelization factors commensurate with full-band or even partial-band reception of any of the civil GNSS signals listed in FIG. 3.

Extension to More Complex Channels

[0184] This signal description, or model, is also now more effectively and easily extensible to more complex channels. Two simple extensions include local multipath channels, in which the simple direct path link model shown in FIG. 1 is extended to include effects of multipath due to scattering or specular reflectors local to and traveling with the receiver platform, and unsynchronized reception scenarios in which the clock driving the ADC and LO has a rate and timing phase that is not fully synchronized to universal time coordinates (UTC). In the former case, Eqn (8) can be extended to

[00031]$\begin{array}{cc}{s}_R\left(t_0+t;\right)\approx e^{j.Math..Math.2.Math.\pi \left({\alpha }_{\mathrm{TR}}\left(\right).Math.t+\frac{1}{2}.Math.{\alpha }_{\mathrm{TR}}^{\left(1\right)}\left(\right).Math.t^2\right)}.Math.g_{\mathrm{TR}}\left(\right).Math.\left(h_R\left(t;\right)*s_T\left(\left(1-{\tau }_{\mathrm{TR}}^{\left(1\right)}\left(\right)\right).Math.t-{\tau }_{\mathrm{TR}}\left(\right);\right)\right)& \mathrm{Eqn}.Math..Math.\left(46\right)\end{array}$

where h.sub.R(t; l) is the impulse response of the channel local to the receiver and “*” denotes the convolution operation. The GNSS link l channel model is then adjusted by replacing h.sub.T(t; l) with

[00032]$\begin{array}{cc}{h}_{\mathrm{TR}}\left(t,\right)=\frac{1}{1-{\tau }_{\mathrm{TR}}^{\left(1\right)}\left(\right)}.Math.\int h_R\left(\frac{u}{1-{\tau }_{\mathrm{TR}}^{\left(1\right)}\left(\right)};\right).Math.h_T\left(t-u,\right).Math.\mathrm{du}.Math.H_{\mathrm{TR}}\left(f;\right)=H_R\left(\left(1-{\tau }_{\mathrm{TR}}^{\left(1\right)}\left(\right)\right).Math.f;\right).Math.H_T\left(f;\right),& \mathrm{Eqn}.Math..Math.\left(47\right)\end{array}$

in Eqn (39), where H.sub.T(f; l), H.sub.R(f; l), and H.sub.TR(f; l) are the Fourier transforms of h.sub.T(t; l), h.sub.R(t; l), and h.sub.TR(t; l), respectively.

[0185] In the latter case, assuming fractional clock rate offset of ε.sub.R and timing offset of t.sub.R, such that the LO downconversion frequency is given by f.sub.R=(1+ε.sub.R)f.sub.T and the ADC sampling times are given by t(n.sub.ADC)=t.sub.R+T.sub.ADCn.sub.ADC, where T.sub.ADC=T.sub.sym/(1+ε.sub.R)M.sub.ADC, then the description derived here is revised by replacing φ.sub.TR(l), τ.sub.TR(l), τ.sub.TR.sup.(1)(l), α.sub.TR(l), and α.sub.TR.sup.(1)(l) with

[00033]$\begin{array}{cc}{\varphi }_{\mathrm{TR}}\left(\right)\leftarrow {\varphi }_{\mathrm{TR}}\left(\right)+2.Math.\pi \left(\left({\alpha }_{\mathrm{TR}}\left(\right)-{\mathrm{.Math.}}_R.Math.f_R\right).Math.t_R+\frac{1}{2}.Math.{\alpha }_{\mathrm{TR}}^{\left(1\right)}\left(\right).Math.t_R^2\right),& \mathrm{Eqn}.Math..Math.\left(48\right)\\ {\tau }_{\mathrm{TR}}\left(\right)\leftarrow {\tau }_{\mathrm{TR}}\left(\right)+t_R.Math.{\tau }_{\mathrm{TR}}^{\left(1\right)}\left(\right),& \mathrm{Eqn}.Math..Math.\left(49\right)\\ {\tau }_{\mathrm{TR}}^{\left(1\right)}\left(\right)\leftarrow \frac{{\tau }_{\mathrm{TR}}^{\left(1\right)}\left(\right)+{\mathrm{.Math.}}_R}{1+{\mathrm{.Math.}}_R},& \mathrm{Eqn}.Math..Math.\left(50\right)\\ {\alpha }_{\mathrm{TR}}\left(\right)\leftarrow \frac{{\alpha }_{\mathrm{TR}}\left(\right)+{\alpha }_{\mathrm{TR}}^{\left(1\right)}\left(\right).Math.t_R-{\mathrm{.Math.}}_R.Math.f_R}{1+{\mathrm{.Math.}}_R},& \mathrm{Eqn}.Math..Math.\left(51\right)\\ {\alpha }_{\mathrm{TR}}^{\left(1\right)}\left(\right)\leftarrow \frac{{\alpha }_{\mathrm{TR}}^{\left(1\right)}\left(\right)}{1+{\mathrm{.Math.}}_R},& \mathrm{Eqn}.Math..Math.\left(52\right)\end{array}$

in all channel descriptions given after the ADC operation.

Extension to Multifeed Reception

[0186] The single-feed reception can be extended simply to multifeed reception shown in FIG. 16, in which the GNSS signals are received over multiple spatial or polarization diverse antennas. In the exemplary direct-conversion system shown in FIG. 16, the GNSS signals are received over M.sub.food feeds using a set of spatial and/or polarization diverse antennas (90a-90c), in some embodiments coupled from a larger set of antennas to M.sub.feed feeds using a beamforming network (BFN) (94). The M.sub.feed feeds are then coherently converted to complex-baseband using a direct-conversion receiver bank (101), to form M.sub.feed×1 signal vector x.sub.R(t−τ.sub.R)=[x.sub.R(t−τ.sub.R; m.sub.feed)].sub.m.sub.feed.sub.=1.sup.M.sup.feed. The feed vector x.sub.R(t−τ.sub.R) is then sampled by a bank of dual-ADC's 143, at sampling rate M.sub.ADCf.sub.sym determined by a common clock (117) where f.sub.sym is the baseband symbol rate of the MOS-DSSS civil GNSS signals transmitted to the receiver, and M.sub.ADC is a positive integer, and channelized into an M.sub.chn×1 signal vector x.sub.chn(n.sub.sym; m.sub.feed) using a bank of symbol-rate synchronous 1:M.sub.chn vector channelizers (160), for example, using a 1:M.sub.ADC S/P operation and an M.sub.chn×M.sub.ADC channelization matrix T.sub.chn. The channelizer output signals) {x.sub.chn(n.sub.sym; m.sub.feed)}.sub.m.sub.feed.sub.=1.sup.M.sup.feed are then combined to form an M.sub.DoF×1 complex vector x.sub.DoF(n.sub.sym) (159), for example, by stacking them over receiver feeds as shown in Eqn (3), where M.sub.DoF=M.sub.feedM.sub.chn is the total number of receiver DoF's used by the invention.

[0187] Assuming coherent downconversion of the received signals, the complex M.sub.feed×1 signal generated at the input to the dual-ADC band is modeled as

[00034]$x_R\left(t-{\tau }_R\right)=i_R\left(t-{\tau }_R\right)+\underset{=1}{\overset{L}{.Math.}}.Math.s_R\left(t-{\tau }_R;\right),$

where i.sub.R(t) comprises the M.sub.feed×1 vector of background noise and CCI present in the GNSS band, and s.sub.R(t−τ.sub.R; l) is M.sub.feed×1 vector of GNSS signals received from transmitted l. In the absence of nonlocal multipath, and over narrow receiver bandwidth, s.sub.R(t; l) is modeled by

[00035]$\begin{array}{cc}{s}_R\left(t_0+t;\right)\approx a_{\mathrm{TR}}\left(t;\right).Math.e^{j.Math..Math.2.Math.\pi \left({\alpha }_{\mathrm{TR}}\left(\right).Math.t+\frac{1}{2}.Math.{\alpha }_{\mathrm{TR}}^{\left(1\right)}\left(\right).Math.t^2\right)}.Math.g_{\mathrm{TR}}\left(\right).Math.s_T\left(\left(1-{\tau }_{\mathrm{TR}}^{\left(1\right)}\left(\right)\right).Math.t-{\tau }_{\mathrm{TR}}\left(\right);\right),& \mathrm{Eqn}.Math..Math.\left(53\right)\end{array}$

where a.sub.TR(t; l) is the M.sub.feed×1 time-varying spatial signature of the signal, which is also assumed to adhere to a first-order DOA blur model due to movement of the receiver over the reception interval.

[0188] Assuming the antennas have nonzero gain along right-hand and left-hand circular polarizations, and in absence of local scattering multipath, a.sub.TR(t; l) can be modeled as

a.sub.TR(t; l)=A.sub.R(Ψ.sub.R(t)u.sub.TR(t; l))ρ.sub.Tr(t;l),   Eqn (54)

where {A.sub.R(u)}.sub.∥u∥.sub.2.sub.=1={[a.sub.R(u)|.sub.RHCPa.sub.R(u)|.sub.LHCP]}.sub.∥u∥.sub.2.sub.=1 is the array manifold of M.sub.feed×2 complex gains along the right-hand and left-hand circular polarizations at the BFN output, parameterized with respect to 3×1 unit-norm local direction-of-arrival (DOA) vector u, and where ρ.sub.TR(t; l) is a (generally time-varying) 2×1 unit-norm polarization gain vector. The array manifold can include adjustments to account for multipath local to the receiver platform, mutual coupling between antenna elements, and direction-independent complex gains in any distribution system coupling the array to the BFN and/or the receiver bank.

[0189] As discussed for the scenario shown in FIG. 9 and summarized in FIG. 11, the local signal DOA's adhere's closely to a first-order model over time intervals on the order of 10 seconds. For an antenna array with close element spacing, the time-varying spatial signature should therefore also adhere to a first-order model, i.e., a.sub.TR(t; l)≈a.sub.TR(t.sub.0; l)+a.sub.TR.sup.(1)(t.sub.0; l)(t−t.sub.0) for |t−t.sub.0|≦10 s, where a.sub.TR.sup.(1)(t; l) da.sub.TR(t; l)/dt is the first derivative of a.sub.TR(t; l). Using this approximation, then the GNSS signal l MIMO link matrix is given by

A.sub.DoF(n.sub.sym;l)=(a.sub.TR(t.sub.0;l)+a.sub.TR.sup.(1)(t.sub.0;l)T.sub.symn.sub.sym)A.sub.chn(n.sub.sym;l)+a.sub.TR.sup.(1)(t.sub.0;l)T.sub.symA.sub.chn.sup.(1)(n.sub.sym;l),   Eqn (55)

where

[00036]$A_{\mathrm{chn}}\left(n_{\mathrm{sym}};\right)=T_{\mathrm{chn}}.Math.A_{\mathrm{sym}}\left(n_{\mathrm{sym}};\right).Math..Math.\mathrm{and}.Math..Math.A_{\mathrm{chn}}^{\left(1\right)}\left(n_{\mathrm{sym}};\right)=T_{\mathrm{chn}}.Math.\mathrm{diag}.Math.\left\{\frac{m_{\mathrm{ADC}}}{M_{\mathrm{ADC}}}\right\}.Math.A_{\mathrm{sym}}\left(n_{\mathrm{sym}};\right),$

and where A.sub.sym(n.sub.sym; l) is given by Eqn (29)-Eqn (31). Equation Eqn (55) also admits low-order regression model Eqn (40) for A.sub.DoF(n.sub.sym; l), allowing s.sub.DoF(n.sub.sym; l) to be modeled by Eqn (42)-Eqn (44). This receiver model also extends to dispersive channels in which the feeds are subject to local multipath, including local multipath that is substantively different on each feed, using channel modifications given in Eqn (46)-Eqn (47); and to receivers with clocks that are not synchronized to UTC, using channel modifications given in Eqn (48)-Eqn (52).

[0190] Note that the structure of the end-to-end multifeed MIMO network model is identical to the structure of the single-feed MIMO network model given in Eqn (33) and Eqn (45). Any differences lie only in the internal structure of the MIMO link matrix A.sub.DoF(n.sub.sym; l) and A′.sub.DoF(l), and the structure of the co-channel interference added to those signals. Therefore, any algorithm that blindly exploits the network model given in Eqn (33) or Eqn (45) can be applied to either a single-feed or a multifeed receiver with no change in structure, albeit with differences in implementation complexity if M.sub.DoF varies as a function of M.sub.feed, and with differences in performance, e.g., if the additional spatial/polarization of the received signals allows the processor to better separate GNSS signals or remove CCI.

[0191] To overcome this complexity issue, the baseline channelizer method can be adjusted to provide the same number of DoF's for the multifeed receiver as the single-feed receiver, by reducing M.sub.chn or the ADC sampling rate by 1/M.sub.feed, such that M.sub.DoF remains constant as M.sub.feed is increased. This is easily accomplished using the thinned-FFT channelizer shown in FIG. 14, by simply reducing the number of output bins used by the channelizer. If the channelizer bins remain densely spaced, this approach has the added benefit of reducing the effect of TOA blur on the output signal. Alternately, if the channelizer bins are sparsely spaced, this approach can reduce complexity of the channelization operation itself. However, it should be noted that this processor sacrifices despreading gain that could improve the SINR of the output symbol stream.

[0192] Importantly, while the TOA and FOA of a GNSS transmitter can be easily spoofed in a covert or “aligned” spoofing scenario, the DOA (and, to a lesser degree, the polarization) of that transmitter cannot be easily spoofed. In addition, the multifeed receiver can null any CCI impinging on the array, if the array has sufficient degrees of freedom to separate that CCI from the GNSS signals.

Blind Civil GNSS Despreading Method

Adaptive Despreader Structures

[0193] The single-feed and multifeed adaptive despreader structures are shown in FIG. 17 and FIG. 18. Both structures comprise the following. [0194] A receiver front-end, which downconverts, digitizes, and channelizes signals received over one or more antennas, providing a M.sub.DoF×1 vector signal x.sub.DoF(n.sub.sym) with sample rate f.sub.sym equal to the GNSS signal symbol rates. [0195] A blind despreading stage, which directly exploits the general MIMO network structure given in Equation Eqn (33) or Eqn (45) to detect, despread, and demodulate received GNSS signals without explicitly using or searching over ranging codes or (for multifeed receivers) the array manifold of the receiver, using only that MIMO network structure and properties of the GNSS symbol streams listed in FIG. 3. Among other advantages, this stage can operate under potentially heavy CCI, including wideband CCI in multifeed receivers, and can greatly exceed the performance of conventional correlative matched-filter despreaders. This stage also estimates the number of transmitted GNSS signals L.sub.T in the receiver's field of view, and estimates key geo-observables of those signals, e.g., the fractional fine FOA and symbol-normalized DFOA (referred to here as the FOA vectors {{tilde over (α)}.sub.TR(l)}.sub.l=1.sup.L.sup.T), initial coarse TOA symbol offsets {n.sub.TR(l)}.sub.l=1.sup.L.sup.T, and in some embodiments MIMO signature matrices {A′.sub.DoF(l)}.sub.l=1.sup.L.sup.T for the detected GNSS signals, as part of the blind despreading procedure. [0196] A optional copy-aided analysis stage, which determines remaining geo-observables needed for a full PNT solution, using copy-aided parameter estimation methods that exploit results of the blind despreading stage and detailed structure of the MIMO link matrices given in Eqn (28)-Eqn (31) for the single-feed system shown in FIG. 17, and the additional structure given in Eqn (54)-Eqn (55) for the multifeed system shown in FIG. 18. The copy-aided analysis stage does require search over code timing to complete the PNT, as well as the array manifold if DOA is accessible and desired, e.g., to unambiguously identify spoofers. However this search is simplified by geo-observables already estimated in the blind despreader, and does not require a search over code PRN if that is accessible within the demodulated navigation signals. Moreover, the copy-aided estimator outperforms estimators obtained using correlative methods, especially in presence of strong CCI and GNSS signals, e.g., co-channel beacons, because those prior statistics are both robust to CCI and cross-interference between other GNSS signals. This performance improvement can be substantive if co-channel beacons are received at strong signal-to-noise ratio (SNR), as in this case the performance of correlative methods will be gated by cross-interference between beacons rather than background noise and interference.

[0197] The PNT solution can also be used to provide clock offset estimates to synchronize the receiver(s) to UTC, possibly on a per-feed basis if the receivers in FIG. 18 are driven by independent clocks (not recommended). In this case, the corrections can be applied directly to the receiver(s), or in digital post-processing after the ADC.

[0198] The despreading structures can be adapted on either a continuous basis, in which geo-observables are updated rapidly over time, or on a batch processing basis, in which a block of N.sub.sym channelized data symbols {x.sub.DoF(n.sub.sym)}.sub.n.sub.sym.sub.=0.sup.N.sup.sym.sup.−1 are computed and passed to a DSP processing element that detects the GNSS signals within that data block. The latter approach is especially useful if the invention is being used to develop resilient PNT analytics to aid a primary navigation system, e.g., to assess quality and availability of new GNSS transmissions, or to detect or confirm spoofing transmissions on a periodic basis. The batch adaptive despreading algorithms are described in more detail below.

Batch Adaptive Despreading Procedure for General Civil GNSS Signals

[0199] An exemplary adaptation procedure generally applicable to batch processing of any of the MOS-DSSS civil GNSS signals described in FIG. 3, and that improves any set of the efficiency and speed of the signal processing, is shown in FIG. 19 and described in more specific detail in the detailed description for that drawing. The algorithm is implemented on a “batch processing” basis, by collecting channelizer output data over N.sub.sym data symbols, and detecting and extracting the civil GNSS navigation streams directly from that data set. The algorithm is designed to perform this processing from a “cold start,” (218) i.e., with no prior information about the signals contained within that data set. However, if the prior FOA's and DFOA's of the signals can be provided to the processor (222), or if baseband symbols or navigation data for the signals can be provided (222), then the procedure can be started at an intermediate point in the processing.

[0200] The first 8 steps of FIG. 20 are those operations performed in the blind despreading stage. The last two steps of FIG. 20 are the copy-aided analysis stage. Note that neither ranging code nor array calibration data is required until the copy-aided analysis stages; in particular, the navigation data can be demodulated without ever searching over that data. This is a marked departure from both conventional correlative or despreading methods in the prior art, which require identification and search over the ranging code before the navigation signal can be demodulated, and “codeless despreading” algorithms that cannot demodulate that data.

[0201] This procedure enables a great deal of refinement and accurate discrimination to more closely constrain and limit the processing necessary to accurately interpret the signal's content, before the copy-aided analysis phase begins. In some use scenarios, the blind despreading stage can in fact obviate the copy-aided analysis phase, e.g., if the invention is developing resilient PNT analytics to aid a primary navigation system, or it can be used to substantively thin the number of transmissions that must be analyzed. This procedure can thereby reduce the processing complexity and considerable feedback lag, enabling quicker, more effective signal discrimination without requiring the full processing and analysis of the signal be completed first (or even together).

Batch Adaptive Despreadinq Procedure Implementation for GPS C/A Signals

[0202] An exemplary adaptation procedure applicable to the GPS C/A signal, and that improves any set of the efficiency and speed of the signal processing is shown in FIG. 20 and described in more specific detail in the detailed description for that drawing. The algorithm is implemented on a “batch processing” basis, by collecting channelizer output data over N.sub.sym data symbols, and detecting and extracting the GPS navigation streams directly from that data set. The algorithm is designed to perform this processing from a “cold start,” i.e., with no prior information about the signals contained within that data set. However, if the prior FOA's of the signals (and in particular FOA's derived from the FOA vectors for those signals) are known, then the procedure can be started at an intermediate point in the processing.

[0203] The first 8 steps of FIG. 20 are those operations performed in the blind despreading stage. The last two steps of FIG. 20 are the copy-aided analysis stage. Note that neither ranging code nor array calibration data is required until the copy-aided analysis stages; in particular, the navigation data can be demodulated without ever searching over that data. This is a marked departure from both conventional correlative or despreading methods in the prior art, which require identification and search over the ranging code before the navigation signal can be demodulated, and “codeless despreading” algorithms that cannot demodulate that data.

[0204] This procedure enables a great deal of refinement and accurate discrimination to more closely constrain and limit the processing necessary to accurately interpret the signal's content, before the copy-aided analysis phase begins. In some use scenarios, the blind despreading stage can in fact obviate the copy-aided analysis phase, e.g., if the invention is developing resilient PNT analytics to aid a primary navigation system, or it can be used to substantively thin the number of transmissions that must be analyzed. This procedure can thereby reduce the processing complexity and considerable feedback lag, enabling quicker, more effective signal discrimination without requiring the full processing and analysis of the signal be completed first (or even together).

Data Whitening Operation

[0205] The batch processing procedures shown in FIG. 19 and FIG. 20 perform a preliminary data whitening operation (220) to simplify subsequent despreading and linear algebraic adaptation algorithms. An exemplary whitening procedure is the QR decomposition (QRD), known to those skilled in that field. Forming N.sub.sym×M.sub.DoF windowed data matrix

[00037]$\begin{array}{cc}{X}_{\mathrm{DoF}}.Math.\stackrel{\Delta }{=}.Math.\left(\begin{array}{c}\sqrt{{\omega }_{\mathrm{DoF}}\left(0\right)}.Math.x_{\mathrm{DoF}}^H\left(0\right)\\ \mathrm{.Math.}\\ \sqrt{{\omega }_{\mathrm{DoF}}\left(N_{\mathrm{sym}}-1\right)}.Math.x_{\mathrm{DoF}}^H\left(N_{\mathrm{sym}}-1\right)\end{array}\right),& \mathrm{Eqn}.Math..Math.\left(56\right)\end{array}$

where {ω.sub.DoF(n.sub.sym)}.sub.n.sub.sym.sub.=0.sup.N.sup.sym.sup.−1 is a real data window satisfying

[00038]$\underset{n_{\mathrm{sym}}=0}{\overset{N_{\mathrm{sym}}-1}{.Math.}}.Math.{\omega }_{\mathrm{DoF}}\left(n_{\mathrm{sym}}\right)=1.Math..Math.\mathrm{and}.Math..Math.{(.Math.)}^H$

denotes the conjugate-transpose (Hermitian) operation, then the QRD of X.sub.DoF, denoted {Q.sub.DoF, R.sub.DoF}=QRD(X.sub.DoF), solves

[00039]$\begin{array}{cc}{R}_{\mathrm{DoF}}=\mathrm{chol}.Math.\left\{X_{\mathrm{DoF}}^H.Math.X_{\mathrm{DoF}}\right\},& \mathrm{Eqn}.Math..Math.\left(57\right)\\ Q_{\mathrm{DoF}}=X_{\mathrm{DoF}}.Math.R_{\mathrm{DoF}}^{-1}=\left(\begin{array}{c}{q}_{\mathrm{DoF}}^H\left(0\right)\\ \mathrm{.Math.}\\ q_{\mathrm{DoF}}^H\left(N_{\mathrm{sym}}-1\right)\end{array}\right),& \mathrm{Eqn}.Math..Math.\left(58\right)\end{array}$

where chol(•) is the Cholesky factorization operation, such that X.sub.DoF=Q.sub.DoFR.sub.DoF and Q.sub.DoF.sup.HQ.sub.DoF=I.sub.M.sub.DOF, where I.sub.M is the M×M identity matrix. The data window is norminally rectangular; however, other windows are recommended if the FOA offset between GNSS signals is small, e.g., for fixed ground beacons or aligned spoofers.

Blind GPS C/A Despreading Algorithms

[0206] The blind GPS C/A despreading flow diagram shown in FIG. 20 comprises the following key processing stages: [0207] An emitter detection operation, which detects GNSS emitters in the environment using a phase-SCORE algorithm that exploits the high single-lag correlation of the GNSS symbol streams, along with additional analysis operations to associate phase-SCORE solutions with specific emitters. [0208] A FOA refinement operation, which refines the fine fractional FOA estimate provided by phase-SCORE, and prunes false alarms, using a frequency doubler spectrum that exploits the BPSK structure of the GNSS symbol streams. [0209] A despreader refinement operation, which refines the phase-SCORE based emitter solution, using a conjugate-SCORE algorithm that also exploits the BPSK structure of the GNSS symbol streams. [0210] A NAV stream demodulation operation, which identifies the 50 sps NAV stream timing phase within each 1 ksps GNSS symbol stream, and demodulates each NAV stream, based on the 1:20 interpolated BPSK structure of the NAV stream.

[0211] These operations are described in more detail in the subsections below.

Phase-SCORE Detection Implementation

Phase-SCORE Eigenequation

[0212] The phase-SCORE implementation exploits the autocorrelation function R.sub.d.sub.T.sub.d.sub.T(m.sub.lag)=d.sub.T(n.sub.sym)d*.sub.T(n.sub.sym−m.sub.lag) induced by the 1:M.sub.NAV interpolation of the transmitted (BPSK) GPC C/A symbol stream, i.e., d.sub.T(n.sub.sym)=b.sub.T(└n.sub.sym/M.sub.NAV┘), given by

[00040]$\begin{array}{cc}{R}_{d_T.Math.d_T}\left(m_{\mathrm{lag}}\right)=\{\begin{array}{cc}1-\frac{.Math.m_{\mathrm{lag}}.Math.}{M_{\mathrm{NAV}}},& -M_{\mathrm{NAV}}<m_{\mathrm{lag}}<M_{\mathrm{NAV}}\\ 0,& \mathrm{otherwise}\end{array}& \mathrm{Eqn}.Math..Math.\left(59\right)\end{array}$

for general 1:M.sub.NAV interpolation, where • denotes infinite time-averaging over n.sub.sym and (•)* denotes the conjugation operation (unused here for BPSK symbol sequences). The phase-SCORE implementation exploits this property, by solving the generalized eigenequation

λ.sub.PSC(m){circumflex over (R)}.sub.x.sub.DoF.sub.x.sub.DoFw.sub.PSC(m)={circumflex over (R)}.sub.x.sub.DoF.sub.x.sub.DoF(m.sub.lag)w.sub.PSC(m), |λ.sub.PSC(m)|≧|λ.sub.PSC(m+1)|,   Eqn (60)

where {circumflex over (R)}.sub.x.sub.DoF.sub.x.sub.DoF (m.sub.lag) is the lag-m.sub.lag sample data autocorrelation matrix (ACM) computed over data set

[00041]$\begin{array}{cc}{\left\{x_{\mathrm{DoF}}\left(n_{\mathrm{sym}}\right)\right\}}_{n_{\mathrm{sym}}=0}^{N_{\mathrm{sym}}-1},.Math.{\hat{R}}_{x_{\mathrm{DoF}}.Math.x_{\mathrm{DoF}}}\left(m_{\mathrm{lag}}\right).Math.\stackrel{\Delta }{=}.Math.\underset{n_{\mathrm{sym}}=m_{\mathrm{lag}}}{\overset{N_{\mathrm{sym}}-1}{.Math.}}.Math.\sqrt{{\omega }_{\mathrm{DoF}}\left(n_{\mathrm{sym}}\right).Math.{\omega }_{\mathrm{DoF}}\left(n_{\mathrm{sym}}-m_{\mathrm{lag}}\right)}.Math.x_{\mathrm{DoF}}\left(n_{\mathrm{sym}}\right).Math.x_{\mathrm{DoF}}^H\left(n_{\mathrm{sym}}-m_{\mathrm{lag}}\right)=X_{\mathrm{DoF}}^H.Math.S\left(m_{\mathrm{lag}}\right).Math.X_{\mathrm{DoF}},.Math.S\left(m_{\mathrm{lag}}\right)=\left(\begin{array}{cc}0& 0\\ I_{N_{\mathrm{sym}}-m_{\mathrm{lag}}}& 0\end{array}\right),& \mathrm{Eqn}.Math..Math.\left(61\right)\end{array}$

and where zero-lag ACM {circumflex over (R)}.sub.x.sub.DoF.sub.x.sub.DoF {circumflex over (R)}.sub.x.sub.DoF.sub.x.sub.DoF(0) and X.sub.DoF is given by Eqn (56).

[0213] Using the QRD described in Eqn (56)-Eqn (58), the phase-SCORE eigenequation given in Eqn (60) is compactly expressed as

[00042]$\begin{array}{cc}{\lambda }_{\mathrm{PSC}}\left(m\right).Math.u_{\mathrm{PSC}}\left(m\right)={\hat{R}}_{q_{\mathrm{DoF}}.Math.q_{\mathrm{DoF}}}\left(m_{\mathrm{lag}}\right).Math.u_{\mathrm{PSC}}\left(m\right),.Math..Math.{\lambda }_{\mathrm{PSC}}\left(m\right).Math.\ge .Math.{\lambda }_{\mathrm{PSC}}\left(m+1\right).Math.,& \mathrm{Eqn}.Math..Math.\left(62\right)\\ \begin{array}{c}{\hat{R}}_{q_{\mathrm{DoF}}.Math.q_{\mathrm{DoF}}}\left(m_{\mathrm{lag}}\right)=.Math.\underset{n_{\mathrm{sym}}=1}{\overset{N_{\mathrm{sym}}-1}{.Math.}}.Math.q_{\mathrm{DoF}}\left(n_{\mathrm{sym}}\right).Math.q_{\mathrm{DoF}}^H\left(n_{\mathrm{sym}}-m_{\mathrm{lag}}\right)\\ =.Math.Q_{\mathrm{DoF}}^H.Math.S\left(m_{\mathrm{lag}}\right).Math.Q_{\mathrm{DoF}},\end{array}& \mathrm{Eqn}.Math..Math.\left(63\right)\end{array}$

where u.sub.PSC(m)=R.sub.DoFw.sub.PSC(m)/∥R.sub.DoFw.sub.PSC(m)∥.sub.2 is the unit-norm whitened phase-SCORE eigenvector for mode m, which forms mode m windowed despread output signal z.sub.PSC(n.sub.sym; m)=u.sub.PSC.sup.H(m)q.sub.DoF(n.sub.sym)=√{square root over (ω.sub.DoF(n.sub.sym))}y.sub.PSC(n.sub.sym; m), and where y.sub.PSC(n.sub.sym; m)=w.sub.PSC.sup.H(m)x.sub.DoF(n.sub.sym). The whitening operation efficiently implements many additional operations performed in the system, reducing computational complexity and load.

[0214] In the presence of TOA blur, and under the same conditions described above, the phase-SCORE eigenequation provides Q′.sub.net signal-capture solutions, separated into L.sub.T sets of eigenequations for each GNSS signal. The signal-capture set for GNSS signal l comprises Q′.sub.sym(l) solutions in which the eigenvalue phase is approximately proportional to the fine fractional FOA {tilde over (α)}.sub.TR(l), and the eigenvector nulls the other GNSS symbol streams {d′.sub.R(n.sub.sym; l′)}.sub.l′≠l={d.sub.R(n.sub.sym; l′)g.sub.blur(n.sub.sym; l′)}.sub.l′≠l, suppresses the received non-GNSS interference, and extracts a linear combination of the elements of d′.sub.R(n.sub.sym; l)=d.sub.R(n.sub.sym; l)g.sub.blur(n.sub.sym; l). This allows phase-SCORE modes to be associated with specific emitters, based on the phase of the phase-SCORE eigenvalues. In practice, phase-SCORE should detect 1+Q.sub.blur(l) modes per GNSS signal, associated with the single-lag symbol stream.

[0215] In the presence of FOA blur due to nonzero DFOA, the lag correlation statistic exploited by the eigenequation is degraded. However, if |{tilde over (α)}.sub.TR.sup.(1)m.sub.lag|N.sub.sym<<1, this degradation will be small. Examination of FIG. 10B for the exemplary scenario shown in FIG. 9 shows that |{tilde over (α)}.sub.TR.sup.(1)|N.sub.sym≦0.03 over reception intervals on the order of 10 s (N.sub.sym=10,000), and for m.sub.lag=1, and this condition holds.

[0216] The autocorrelation function given in Equation Eqn (59) is a well known property of the GPS C/A signal, and has been explicitly exploited to design algorithms for adaptation of spatial antijam arrays using the closely-related cross-SCORE algorithm. Yet, while as shown in the prior art the cross-SCORE algorithm also possesses an eigensorting property, it fails to sort signals with the same correlation strength. In contrast, the phase-SCORE algorithm can eigensort signals with differing autocorrelation phase, even if they have the same strength. And, if the receiver front-end is spatially diverse, it can provide the same ability to null wideband CCI as the spatial antijam methods in the prior art.

[0217] Note that the GNSS signal capture solutions also automatically reject narrowband CCI (NBCCI), if that NBCCI occupies a small number of output channels. However, in the limiting case where the NBCCI is a CW tone, the phase-SCORE algorithm will also possess a single tone capture solution that extracts that tone from the received data. In that case, additional screening algorithms must used to identify and remove those tones ahead of any subsequent demodulation operations. Also, even in that case, the system should possess sufficient degrees of freedom to remove literally dozens of those tones from the GNSS signal capture solutions.

[0218] With the geo-observables (Frequency, Timing, and Direction of Arrival) distinguished for a received signal, any set of those values can be compared against a set for such which are known to be correct, i.e. associated with an authorized transmitter—just as a phone number, or Caller ID text link, can be compared against a listing of authorized callers—all without reference to the content of the incoming signal, or message. Signals with the wrong metadata can then be dropped from any further processing, improving the efficiency of the receiver. Only those with correct metadata or not-proven-wrong-yet metadata, need to be further processed to examine the content and structure of the message. This is akin to hearing a shout across a field from a direction (or time) where (or when) no friendly speaker can be calling; or receiving a phone call, text, or email purporting to be from the IRS but originating in Ghana. You don't need to understand any such message to know it's a fake.

Phase-SCORE Detection Operations

[0219] The phase-SCORE detection operations comprise the following steps: [0220] A whitened correlation matrix computation operation to compute M.sub.DoF×M.sub.DoF matrix {circumflex over (R)}.sub.q.sub.DoF.sub.q.sub.DoF(1) using Eqn (63). [0221] An eigenvalue computation operation to compute eigenvalues {λ.sub.PSC(m)}.sub.m=1.sup.M.sup.DoF solving Eqn (62), using mature methods such as the QR algorithm. [0222] A preliminary mode selection operation to identify M.sub.PSC candidate signal-capture solutions based on {|λ.sub.PSC(m)|}.sub.m=1.sup.M.sup.DoF. [0223] An eigenvector computation operation to compute eigenvectors {u.sub.PSC(m)}.sub.m=1.sup.M.sup.PSC for the M.sub.PSC preliminary candidate signal-capture solutions, using mature methods such as Schur decomposition and back-substitution algorithms. [0224] A FOA vector optimization procedure, to estimate FOA vectors {{circumflex over (α)}.sub.PSC(m)}.sub.m=1.sup.M.sub.PSC for each candidate signal-capture solution, using generalized frequency doubler spectra generated from {z.sub.PSC(n.sub.sym; m)}.sub.m=1.sup.M.sub.PSC, where z.sub.PSC(n.sub.sym; m)=u.sub.PSC.sup.H(m)q.sub.PSC(n.sub.sym). [0225] A mode screening procedure, to select {circumflex over (M)}.sub.PSC finalized GNSS signal-capture modes based on the optimized generalized frequency doubler spectra. In some embodiments additional statistics are computed and used to remove modes associated with non-GNSS signals, e.g., tonal interferers. [0226] A mode association procedure, to associate final GNSS signal-capture modes with specific GNSS signals based on commonality of FOA vectors, estimate the number of GNSS signals {circumflex over (L)}.sub.T detected by the algorithm, and to estimate symbol-normalized FOA vectors

[00043]${\left\{{\tilde{\alpha }}_{\mathrm{TR}}\left(\right)\right\}}_{=1}^{{\hat{L}}_T}$

comprising the fine fractional FOA and symbol-normalized DFOA for each detected GNSS signal.

[0227] In some embodiments, additional processing is performed to further optimize GNSS signal capture weights for each detected GNSS signal.

[0228] Exemplary mode selection, FOA optimization, mode screening, and mode association operations are described below.

Preliminary Mode Selection Operations

[0229] Exemplary methods for selecting M.sub.PSC preliminary signal-capture modes include the following: [0230] Mode cutoff sorting procedures, in M.sub.PSC-strongest eigenvalues computed using Eqn (62)-Eqn (63), where M.sub.PSC is a preset cutoff parameter, and {λ.sub.PSC(m), u.sub.PSC(m)}.sub.m=1.sup.M.sup.PSC are selected as candidate modes and sent on to confirmation stages. [0231] Model-fitting procedures, in which M.sub.PSC is adaptively computed by fitting {|λ.sub.PSC(m)|}.sub.m=1.sup.M.sup.DoF against a signal-present and signal-absent profile. A reliable model-fitting approach based on the minimum description length (MDL) principle has been developed and tested using the phase-SCORE eigenvalue strength model

[00044]$\begin{array}{cc}.Math.{\lambda }_{\mathrm{PSC}}\left(m\right).Math.~\{\begin{array}{cc}s\left(m\right)+b^T\left(m\right).Math.c_{\mathrm{fit}}& m=1,\mathrm{.Math.}.Math.,M_{\mathrm{PSC}}\\ b^T\left(m\right).Math.c_{\mathrm{fit}},& m=M_{\mathrm{PSC}}+1,\mathrm{.Math.}.Math.,M_{\mathrm{DoF}}\end{array},& \mathrm{Eqn}.Math..Math.\left(64\right)\end{array}$

where {b(m)}.sub.m=1.sup.M.sup.DoF is an order-Q.sub.fit deterministic basis function, e.g., the polynomial basis b(m)=[(m−1).sup.q].sub.q=0.sup.Q.sup.fit, and (•).sup.T denotes the matrix transpose operation. In one embodiment, the MDL procedure is implemented in the following manner: [0232] Precompute (M.sub.DoF−Q.sub.fit−1)×M.sub.DoF matrix B.sub.⊥=null(B), where B=[b(1) . . . b(M.sub.DoF)].sup.T is the (Q.sub.fit+1)×M.sub.DoF basis matrix and null(•) is the matrix null operation, and pseudoinverse null-matrix set

[00045]${\left\{B_{\perp }^{†}\left(M_{\mathrm{PSC}}\right)\right\}}_{M_{\mathrm{PSC}}=1}^{M_{\mathrm{PSC}}^{\max }},$ [0233] where B.sub.⊥.sup.\(M.sub.PSC) (B.sub.⊥.sup.T(:,1:M.sub.PSC)B.sub.⊥(:,1:M.sub.PSC)).sup.−1B.sub.⊥.sup.T(:,1:M.sub.PSC) is the M.sub.PSC×(M.sub.DoF−Q.sub.fit−1) pseudoinverse of the first M.sub.PSC columns of B.sub.⊥. [0234] Initialize:

[00046]${\mathrm{.Math.}}_{\perp }=B_{\perp }\left(\begin{array}{c}.Math.{\lambda }_{\mathrm{PSC}}\left(1\right).Math.\\ \mathrm{.Math.}\\ .Math.{\lambda }_{\mathrm{PSC}}\left(M_{\mathrm{DoF}}\right).Math.\end{array}\right),\begin{array}{c}\left(M_{\mathrm{DoF}}-Q_{\mathrm{fit}}-1\right)\times 1\\ \mathrm{signal}.Math.\text{-}.Math.\mathrm{absent}.Math..Math.\mathrm{regression}.Math..Math.\mathrm{statistic}\end{array}$$F_{\mathrm{LS}}\left(0\right)={.Math.{\mathrm{.Math.}}_{\perp }.Math.}_2^2,\begin{array}{c}\mathrm{signal}.Math.\text{-}.Math.\mathrm{absent}.Math..Math.\mathrm{least}.Math.\text{-}.Math.\mathrm{squares}\\ \mathrm{regression}.Math..Math.\mathrm{error}\end{array}$$S_{\mathrm{MDL}}\left(0\right)=\frac{Q_{\mathrm{fit}}}{M_{\mathrm{DoF}}}.Math.\ln \left(M_{\mathrm{DoF}}\right)-\ln \left(F_{\mathrm{LS}}\left(0\right)\right),\mathrm{signal}.Math.\text{-}.Math.\mathrm{absent}.Math..Math.\mathrm{MDL}.Math..Math.\mathrm{metric}$ [0235] Iterate until maximized:

[00047]$F_{\mathrm{LS}}\left(M_{\mathrm{PSC}}^{\prime }\right)=F_{\mathrm{LS}}\left(0\right)-{.Math.B_{\perp }^{†}\left(M_{\mathrm{PSC}}^{\prime }\right).Math.{\mathrm{.Math.}}_{\perp }.Math.}_2^2,\begin{array}{c}{M}_{\mathrm{PSC}}^{\prime }.Math.\text{-}.Math.\mathrm{mode}.Math..Math.\mathrm{least}.Math.\text{-}.Math.\mathrm{squares}\\ \mathrm{regression}.Math..Math.\mathrm{error}\end{array}$$S_{\mathrm{MDL}}\left(M_{\mathrm{PSC}}^{\prime }\right)=\frac{M_{\mathrm{PSC}}^{\prime }+Q_{\mathrm{fit}}}{M_{\mathrm{DoF}}}.Math.\ln \left(M_{\mathrm{DoF}}\right)-\ln \left(F_{\mathrm{LS}}\left(M_{\mathrm{PSC}}^{\prime }\right)\right),.Math.M_{\mathrm{PSC}}^{\prime }.Math.\text{-}.Math.\mathrm{mode}.Math..Math.\mathrm{MDL}.Math..Math.\mathrm{metric}$ [0236] The MDL mode loading estimate is then

[00048]$M_{\mathrm{PSC}}=\arg .Math..Math.\underset{M_{\mathrm{PSC}}^{\prime }\in \left\{0,M_{\mathrm{PSC}}^{\mathrm{ma}.Math..Math.x}\right\}}{\max }.Math.S_{\mathrm{MDL}}\left(M_{\mathrm{PSC}}^{\prime }\right),$ [0237] with M.sub.PSC=0 denoting signal absence.

[0238] The mode cutoff procedure works particularly well in environments where the number of potential candidate modes is small and/or predictable, especially if reliable confirmation methods can be used to prune modes developed under the procedure.

FOA Vector Optimization Operations

[0239] Once M.sub.PSC preliminary candidate phase-SCORE signal capture modes have been identified, and the M.sub.PSC whitened eigenvectors {u.sub.PSC(m)}.sub.m=1.sup.M.sub.PSC associated with those modes have been computed, the FOA vectors of the signal captured by each candidate mode are determined by maximizing the magnitude of the generalized frequency-doubler spectrum, given by

[00049]$\begin{array}{cc}\begin{array}{c}{\hat{R}}_{{\mathrm{yy}}^{*}}\left(\alpha \right).Math.\stackrel{\Delta }{=}.Math..Math.\underset{n_{\mathrm{sym}}=0}{\overset{N_{\mathrm{sym}}-1}{.Math.}}.Math.{\omega }_{\mathrm{DoF}}\left(n_{\mathrm{sym}}\right).Math.y^2\left(n_{\mathrm{sym}}\right).Math.\exp \left(-j.Math..Math.4.Math.\pi .Math..Math.g_{\mathrm{FOA}}^T\left(n_{\mathrm{sym}}\right).Math.\alpha \right)\\ =.Math.\underset{n_{\mathrm{sym}}=0}{\overset{N_{\mathrm{sym}}-1}{.Math.}}.Math.z^2\left(n_{\mathrm{sym}}\right).Math.\exp \left(-j.Math..Math.4.Math.\pi .Math..Math.g_{\mathrm{FOA}}^T\left(n_{\mathrm{sym}}\right).Math.\alpha \right),\end{array}& \mathrm{Eqn}.Math..Math.\left(65\right)\end{array}$

for general despreader output signal y(n.sub.sym) w.sup.Hx.sub.DoF (n.sub.sym) and associated whitened despreader output signal z(n.sub.sym)=u.sup.Hq.sub.DoF(n.sub.sym), where u=R.sub.DoFw, and where g.sub.FOA(n.sub.sym)=[g.sub.FOA.sup.(q.sup.FOA.sup.)(n.sub.sym)].sub.q.sub.FOA.sup.Q.sup.FOA is a Q.sub.FOA×1 deterministic FOA basis function, and α=[α(q.sub.FOA)].sub.q.sub.FOA.sup.Q.sup.FOA is a Q.sub.FOA×1 unknown real coefficient vector, referred to here as the FOA vector. If

[00050]$g_{\mathrm{FOA}}^{\left(q_{\mathrm{FOA}}\right)}\left(n_{\mathrm{sym}}\right)=\frac{1}{q_{\mathrm{FOA}}!}.Math.n_{\mathrm{sym}}^{q_{\mathrm{FOA}}},$

referred to here as the polynomial basis, then α(1)={tilde over (α)}.sub.TR (and is the fine fractional FOA), and α(2)={tilde over (α)}.sub.TR.sup.(1) (and is the symbol-normalized first-order DFOA); however, Eqn (65) can be used to estimate the signal FOA and DFOA's of arbitrarily high-order. Eqn (65) can also be defined using other FOA bases, if warranted based on characteristics of the transmission channel. If Q.sub.FOA=1 and g.sub.FOA.sup.(1)(n.sub.sym)=n.sub.sym, then Eqn (65) reduces to the conventional frequency-doubler spectrum described in the prior art, and referred to as the zero-lag conjugate-cyclic autocorrelation function {circumflex over (R)}.sub.yy*.sup.2α described in the prior art as a well-known means for detection and “open-loop” carrier recovery of BPSK and BPSK-DSSS signals.

[0240] The generalized frequency-doubler spectrum can be derived as a maximum-likelihood estimator, by modeling y(n.sub.sym) as

y(n.sub.sym)=ε(n.sub.sym)+d(n.sub.sym)e.sup.jφ(n.sup.sym)   Eq. 66

where ε(n.sub.sym) is an independent, identically-distributed (i.i.d.) circularly-symmetric complex-Gaussian (CSCG) noise sequence with mean zero and unknown nonnegative variance σ.sup.2, d(n.sub.sym) is an unknown real signal sequence, and φ(n.sub.sym)=φ+2πg.sub.FOA.sup.T(n.sub.sym)α is a time-varying phase sequence, and where φ is an unknown real signal phase. The joint ML estimate of α, φ, σ.sup.2, and d(n.sub.sym) under these assumptions is then given by

[00051]$\begin{array}{cc}{\hat{\alpha }}_{\mathrm{ML}}=\arg .Math..Math.\underset{\alpha \in {ℝ}^{Q_{\mathrm{FOA}}}}{\max }.Math..Math.{\hat{R}}_{\mathrm{yy}*}\left(\alpha \right).Math.,& \mathrm{Eqn}.Math..Math.\left(67\right)\\ {\hat{\varphi }}_{\mathrm{ML}}=\mathrm{\pi \delta }+\frac{1}{2}.Math.\arg .Math.{\hat{R}}_{\mathrm{yy}*}\left({\hat{\alpha }}_{\mathrm{ML}}\right),\delta =0,1,& \mathrm{Eqn}.Math..Math.\left(68\right)\\ {\hat{\sigma }}_{\mathrm{ML}}^2=\frac{1}{2}.Math.\left({\hat{R}}_{\mathrm{yy}}-.Math.{\hat{R}}_{\mathrm{yy}*}\left({\hat{\alpha }}_{\mathrm{ML}}\right).Math.\right),& \mathrm{Eqn}.Math..Math.\left(69\right)\\ {\hat{d}}_{\mathrm{ML}}\left(n_{\mathrm{sym}}\right)=\mathrm{Re}.Math.\left\{y\left(n_{\mathrm{sym}}\right).Math.\exp \left(-j\left({\hat{\varphi }}_{\mathrm{ML}}+2.Math.\pi .Math..Math.g_{\mathrm{FOA}}^T\left(n_{\mathrm{sym}}\right).Math.{\hat{\alpha }}_{\mathrm{ML}}\right)\right)\right\}.& \mathrm{Eqn}.Math..Math.\left(70\right)\end{array}$

[0241] Note that {circumflex over (φ)}.sub.ML possesses an ambiguity of 0 or π, resulting in a sign ambiguity of ±1 in {circumflex over (d)}.sub.ML(n.sub.sym). Similarly, the polynomial basis induces an ambiguity of 0 or ±½ in {circumflex over (α)}.sub.ML(1), since {circumflex over (R)}.sub.yy*(α′)={circumflex over (R)}.sub.yy*(α) if g.sub.FOA.sup.(1)(n.sub.sym)=n.sub.sym and

[00052]${\alpha }^{\prime }\left(1\right)=\alpha \left(1\right)+\frac{1}{2},$

resulting in a further frequency modulation ambiguity of 1 or (−1).sup.n.sup.sym in {circumflex over (d)}.sub.ML(n.sub.sym). The frequency modulation ambiguity can be resolved for each mode from the phase of the phase-SCORE eigenvalue {circumflex over (α)}.sub.PSC(m), which is roughly proportional to the true FOA. Additional information is needed to resolve the sign ambiguity, e.g., by identifying the TLM preamble in the GPS navigation message stream.

[0242] In one embodiment, the maximal value of Eqn (65) is robustly found for the polynomial basis and Q.sub.FOA=2 using the following procedure: [0243] Compute the generalized discrete frequency doubler spectrum given by

[00053]$\begin{array}{cc}S\left(k,\right)={.Math.{\mathrm{FFT}}_{K_{\mathrm{FFT}}}\left(z^2\left(n_{\mathrm{sym}}\right).Math.\exp \left(-j.Math..Math.2.Math.{\mathrm{\pi \Delta }}_{\mathrm{DFOA}}.Math.\frac{}{L}.Math.n_{\mathrm{sym}}^2\right)\right).Math.}^2,& \mathrm{Eqn}.Math..Math.\left(71\right)\end{array}$

for l=−L−1, . . . , L+1 and k=0, . . . , K.sub.FFT−1 where FFT.sub.K.sub.FFT(•) is the fast-Fourier transform (FFT) operation with K.sub.FFT≧N.sub.sym output bins (for example, a power-of-two greater than 2N.sub.sym), z(n.sub.sym) is given by Eqn (65), and Δ.sub.FOA covers the range of expected symbol-synchronized DFOA's, e.g., 3.07×10.sup.−6 for the reception scenario shown in FIG. 9, and where k=0, . . . κ.sub.FFT−1 is the output FFT bin index.

[00054]$\mathrm{Set}\left(k_{\mathrm{ma}.Math..Math.x},{}_{\mathrm{ma}.Math..Math.x}\right)=\arg .Math..Math.\underset{k\in \left\{0,\mathrm{.Math.}.Math.,K_{\mathrm{FFT}}-1\right\}.Math..Math.\in \left\{-L,\mathrm{.Math.}.Math.,L\right\}}{\max }.Math.S\left(k,\right),$

and refine each index to sub-bin accuracy using a quadratic-fit algorithm, e.g.,

[00055]$k_{\max }\leftarrow k_{\max }+\frac{1}{2}.Math.\frac{S\left(\left(k_{\max }-1\right).Math.\mathrm{mod}.Math..Math.K_{\mathrm{FFT}},{}_{\max }\right)-S\left(\left(k_{\max }+1\right).Math.\mathrm{mod}.Math..Math.K_{\mathrm{FFT}},{}_{\max }\right)}{\begin{array}{c}S\left(\left(k_{\max }-1\right).Math.\mathrm{mod}.Math..Math.K_{\mathrm{FFT}},{}_{\max }\right)-2.Math.S\left(k_{\max },{}_{\max }\right)+\\ S\left(\left(k_{\max }+1\right).Math.\mathrm{mod}.Math..Math.K_{\mathrm{FFT}},{}_{\max }\right)\end{array}}$${}_{\max }\leftarrow {}_{\max }+\frac{1}{2}.Math.\frac{S\left(k_{\max },{}_{\max }-1\right)-S\left(k_{\max },{}_{\max }+1\right)}{S\left(k_{\max },{}_{\max }-1\right)-2.Math.S\left(k_{\max },{}_{\max }\right)+S\left(k_{\max },{}_{\max }+1\right)}.$