Patent Publication Number: US-11650328-B2

Title: Exploitation of pilot signals for blind resilient detection and geo-observable estimation of navigation signals

Description:
CROSS REFERENCE TO RELATED APPLICATIONS 
     This application is a Continuation of U.S. patent application Ser. No. 16/701,144, filed on Dec. 2, 2019, now U.S. Pat. No. 11,125,888, which claims priority to Provisional Appl. No. 62/773,589, filed Nov. 30, 2018; and Provisional Appl. No. 62/773,605, filed Nov. 30, 2018; and U.S. patent application Ser. No. 16/701,144 is a Continuation-in-Part of U.S. patent application Ser. No. 15/731,417, filed Jun. 5, 2017, now U.S. Pat. No. 10,775,510, which claims priority to Provisional Appl. No. 62/392,623, filed Jun. 6, 2016 and Provisional Appl. No. 62/429,029, filed Dec. 1, 2016; all of which are hereby incorporated by reference in their entireties and all of which this application claims priority under at least 35 U.S.C. 120 and/or any other applicable provision in Title 35 of the United States Code. 
    
    
     FIELD OF THE INVENTION 
     Aspects of the disclosure relate generally to location determination, and more particularly to detecting or determining geo-observables of navigation signals transmitted by satellite vehicles (SV&#39;s) in Global Navigation Satellite Systems (GNSS), airborne non-SV beacons (e.g., pseudolites), and ground-based beacons. 
     BACKGROUND 
     The American Global Positioning System (GPS) L1 coarse acquisition signal (commonly referred to as the “L1 C/A signal,” or more recently as the “L1 legacy signal”), originally developed as a means for aiding acquisition of the military-grade P (later encrypted P(Y)) code, has gained near-ubiquitous use as a means for navigation and carrier/timing synchronization in civil positioning, navigation, and timing (PNT) systems. This is due in part to the GPS network&#39;s worldwide visibility; the publically disclosed and well-known structure of the L1 C/A signal-in-space (SiS); and the ready availability of mature, low-cost hardware for receiving the L1 C/A signal. 
     The success of the GPS L1 C/A signal has resulted in incorporation of “civil” modes in all of the Global Navigation Satellite System (GNSS) signals developed to date, e.g., the GLONASS, Galileo, Beidou, and NAVIC/IRNSS signals. In addition, the need for additional precision to support civil PNT operation in urban canyons and autonomous navigation systems has spurred the development of the wideband civil GPS L5, Galileo E5B, and E6B/C modes, and the Japanese Quasi-Zenith Satellite System (QZSS), which is transmits a version of the civil GPS L5 signal from geosynchronous orbit. Moreover, similar civil signal modes have been proposed as GNSS-like ground and air beacons, such as the Locata&#39;s LocataLite signal, which emulates the GPS C/A signal with a ×10 spreading rate in the 2.4-2.485 GHz ISM band. All of these civil modes possess a common modulation format, referred to herein as modulation-on-symbol direct-sequence spread spectrum (MOS-DSSS), in which a baseband symbol stream, typically (but not always) operating at a 1,000 sample-per-second (1 ksps) symbol rate, is modulated by a higher-rate ranging code that repeats every symbol. 
     Since its inception, the only method used to acquire, demodulate, and determine geo-observables of these navigation signals in PNT applications has been matched-filter despreading methods. These methods operate by correlating the ranging code (or a set of ranging codes that the transmitter may be using) against the received signal and adjusting the carrier and timing offset of that code (referred to herein as geo-observables of the received navigation signal) until it aligns with the code received from the GNSS transmitter. However, matched filter despreading possesses inherent limitations that can degrade its performance in many reception environments. Known limitations include the following:
         Slow cold-start time-to-first-fix (TTFF), even in the most benign environments. In a cold start scenario, the receiver does not know the specific ranging codes used by any of the GNSS satellite vehicles, nor does it know the timing phase or frequency offset of those codes, hence it must search over code, timing, and frequency to synchronize with a single satellite. Cold-start TTFF can be on the order of minutes, and in addition can be quite power consumptive.   High sensitivity to multipath interference, including inadvertent multipath induced by the specular reflectors in the vicinity of the receiver, which can distort the ranging codes and degrade or destroy correlation peaks, and synthetic multipath induced by malicious repeaters. Inadvertent multipath is especially acute in urban and indoor navigation and localization scenarios.   High susceptibility to strong co-channel interference, including near-far interference from local beacons, malicious non-targeted spoofing signals, and malicious jamming signals. The GNSS signals deployed to date only provide a despreading gain of 30-40 dB, based on the chip-rate of the GNSS ranging codes, most of which is needed to raise the GNSS signal above the background noise; hence the signal is easily jammed by many types of interference. An additional 10-13 dB gain can be obtained by exploiting additional structure of the baseband symbol stream. However, this further increases TTFF in cold-start acquisition scenarios.   High susceptibility to targeted or “aligned” covert spoofers that can use knowledge of the ranging code from the satellite and rough location of the receiver platform to duplicate and attempt to supplant signals from those satellites.   Reliance on ancillary anti-spoofing (AS) codes, e.g., the GPS L1 P(Y) code, to overcome jamming and spoofing events. This reliance increases the cost of the receiver, limits AS capability to a subset of military receivers (again, with inherently higher cost), and introduces a host of physical security measures and restrictions to prevent compromise of AS-capable receivers.   Need to maintain a library of all of the GNSS ranging codes, or seeds used to initiate the codes, in order to implement the despreader. This can complicate the receiver architecture, especially for “multi-GNSS receivers” that are designed to exploit civil signals generated by multiple families of satellites.
 
Moreover, there is a perceived (but not real) need for matched-filter despreaders to receive and sample signals over wide bandwidths, constrained to a large integer multiple of the ranging code chip-rate, e.g., 1.023 million chips per second (Mcps) for the GPS C/A code, and 10.23 Mcps for the GPS L5 short code, in order to fully exploit the processing gain of the code, overcome aforementioned channel multipath effects, and (in military GPS receivers) exploit moderate and long codes overlaid on top of the short civil code. This perceived need is further complicated in multi-GNSS radios, which must accommodate wide ranges of chip rates and bandwidths.
       

     SUMMARY 
     The following summarizes some aspects of the present disclosure to provide a basic understanding of the discussed technology. This summary is not an extensive overview of all contemplated features of the disclosure, and is intended neither to identify key or critical elements of all aspects of the disclosure nor to delineate the scope of any or all aspects of the disclosure. Its sole purpose is to present some concepts of one or more aspects of the disclosure in summary form as a prelude to the more detailed description that follows. 
     Disclosed aspects provide for detecting and obtaining geo-observables of navigation signals generated by satellite vehicles (SV&#39;s) in Global Navigation Satellite Systems (GNSS), and generated by ground-based or airborne non-SV beacons, e.g., pseudolites. This includes effecting detection and geo-observable estimation in environments populated with many sources, including malicious sources intending to disrupt or subvert information provided by legitimate sources, for example, non-GNSS jammers, spoofers that emulate GNSS signals, and repeaters that can record and replay environments containing true GNSS signals. Furthermore, this includes effecting said detection and geo-observable estimation in environments subject to potentially severe differences in power levels between such signals, for example, “near-far interference” between remote GNSS signals and nearby ground or air based beacons, jammers, and spoofers. 
     Aspects of the disclosure can achieve these and other goals by exploiting the massive spectral redundancy of navigation signals employing modulation-on-symbol direct-sequence spread spectrum (MOS-DSSS) modulation formats, in which a baseband symbol sequence is spread by a ranging code that repeats every each baseband symbol. Novel and powerful methods for exploiting this signal structure are introduced. For example, aspects described herein can exploit the time, frequency and code diversity of those signals, to implement adaptive linear combining methods and linear-algebraic combiner adaptation algorithms to detect all of the MOS-DSSS navigation signals in the received environment at the maximum signal-to-interference-and-noise ratio (max-SINR) achievable by the combiners; estimate key geo-observables of those detected signals; identify detected malicious sources ahead of or during PNT acquisition and tracking operations; prevent those malicious sources from corrupting or subverting the PNT operation; and optionally alert the receiver to the presence of those sources, and/or geolocate those malicious sources. In aspects employing multifeed receivers wherein multiple spatial and/or polarization diverse antennas are coupled to multiple receivers that are receptive to the MOS-DSSS navigation signals, aspects of the disclosure can additionally exploit the additional spatial and/or polarization diversity of navigation and non-navigation (e.g., jamming) signals to further remove non-navigation signals from the received MOS-DSSS navigation signals, and to remove targeted spoofing signals that closely emulate the time, frequency, and code diversity of navigation signals the spoofers are intending to displace at a victim receiver. 
     Moreover, such aspects can perform all of these operations without knowledge of, a search over, or synchronization to the specific ranging codes employed by each navigation signal. As a consequence, the system described herein can enable faster time-to-first-fix (TTFF) than conventional navigation receivers that require a search over ranging code parameters, and additional robustness to multipath interference (including synthetic multipath induced by repeaters) that can impair that search. Disclosed aspects also provide enhanced compatibility with implementations employing software-defined radio (SDR) architectures, including general-purpose SDR (GPSDR) architectures that are not optimized to GNSS applications, further improving cost and flexibility of the system. 
     Aspects disclosed herein can adapt symbol-rate (e.g., not chip-rate) synchronous reception and channelization structure, using both purpose-built equipment designed specifically for GNSS bands and signals of interest, or general-purpose software-defined receiver (GPSDR) equipment, e.g., Ettus Research Universal Research Software Peripheral (USRP) radios such as the B200mini-series device, which can employ a much wider range of sampling rates and bandwidths, and which employ analog-to-digital converters (ADC) with 8 or more bits per in-phase (I) and quadrature (Q) rail, allowing implementation of sophisticated jammer mitigation techniques and operation in severe multipath environments. This is in contrast to GNSS-specific SDR (GNSS-SDR) equipment, e.g., the Maxim MAX2769 series of Universal GNSS Receivers, which are employ ADC&#39;s with sampling rates preset to integer multiple of the GNSS chip-rate, e.g., 16.238 Msps (16×1.023 Msps), and with ADC precisions of 1-2 bits per I and Q rail, or 3 bits per I rail, thereby substantively worsening the receivers&#39; vulnerability to jamming measures and severe multipath. 
     Aspects disclosed herein can be configured to operate in devices, systems, and methods disclosed in nonprovisional patent application Ser. No. 15/731,417, “Blind Despreading of Civil GNSS Signals for Resilient PNT Applications” (the &#39;417 NPA), incorporated herein by reference. Such aspects can be employed in a method, system, and/or apparatus 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 PNT tracking operations; and for geolocating malicious GNSS signals. Some aspects can be configured for 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. 
     Aspects of the disclosure are further described for exploiting navigation signal modulation formats in which a component of the baseband symbol sequence, referred to herein as a “pilot signal,” is periodically repeated. Novel and powerful methods for exploiting this signal structure is an additional focus of aspects disclosed herein. For example, means are described for exploiting the period or content of those pilot signals to effect “partially blind” signal detection and geo-observable estimation methods that are highly efficient relative to equivalent “fully blind” methods that exploit more general properties of the baseband symbol sequence; and that can perform this detection and geo-observable estimation at much lower received incident powers (RIP&#39;s) than equivalent fully-blind methods. 
     Aspects of the disclosure are further described for exploiting navigation or symbol streams provided by third parties. It is known that fully-blind despreading methods introduce complexity inherent to the goals and needs of the method; specifically, the demodulation (symbol estimation) step. This step inherently requires a despreader output signal-to-noise ratio (SNR) that is high enough to determine the baseband symbols after integration over each symbol, or to demodulate the navigation sequence after integration over each navigation symbol. This introduces a fundamental “acquisition threshold” for any fully blind method, i.e., the received incident power (RIP) or despread signal-to-inference-and-noise ratio (SINR) at which navigation signals can be reliably demodulated by the methods, or for that matter for any conventional matched-filter despreading method. This also places a limit on the accuracy of any geo-observable estimation algorithm that relies on that demodulated navigation data, which will be degraded by errors in that data. 
     However, in many use scenarios, baseband symbol or navigation sequences can be obtained from third-party sources. For example, in some use scenarios, PNT-capable platforms may possess GNSS reception devices that can receive, despread, and demodulate navigation data on their own in non-challenged environments. For example, recent versions of Android provide raw GPS navigation bit-streams along with GNSS measurement results. In other use scenarios, a PNT-capable or noncapable platform may possess communication links to other PNT-capable platforms or network infrastructure capable of receiving, despreading, and demodulating navigation signals. For example CORS (Continuously-Operating Reference Stations) decode navigation data and make it available in a decoded format that can be obtained over the Internet. 
     This third-party data can be used to implement “partially-blind” methods that can achieve the blind signal separation, jammer/spoofer resilience, near-far interference resilience, and multipath robustness without the need for, or sensitivity to, demodulating the received navigation data. Among other advantages, these approaches can develop resilient PNT (RPNT) analytics over integration times that are substantively (in theory, arbitrarily) longer than a single baseband or navigation symbol, allowing accurate signal geo-observable estimation well below the acquisition threshold required for both any fully-blind despreading method. In many use scenarios, this performance can furthermore be provided at substantively lower complexity than fully-blind methods that require demodulation of the navigation sequence. 
     Aspects disclosed herein also allow the use of multi-subband architectures in which symbol-rate synchronous channelizers are employed over multiple subbands within the GNSS signal band, or symbol-rate synchronous channelizer bins covering the GNSS signal bandwidth are organized into subbands within the GNSS signal band. In these aspects, linear combiners, and interference-excising linear-algebraic combiner adaptation algorithms, are employed to develop detection features and signal estimates within each subband metrics, which are then combined across subbands. The resultant method can exploit the full bandwidth of the GNSS signal, with a system complexity that scales linearly with the bandwidth exploited (number of subbands), and can constrain the dimensionality (and therefore complexity) of linear combiner adaptation operations employed in each subband to just that needed to effectively excise interference encountered within that subband. The multi-subband approach also provides a means for naturally excising non-tonal narrowband jamming signals, by selectively and adaptively suppressing subbands containing strong non-GNSS interference. 
     In one aspect, a low-cost single-feed reception system and multi-subband symbol-rate synchronous channelizer is presented that can develop resilient positioning, navigation, and timing (RPNT) analytics for civil modes of Global Navigation Satellite System (GNSS) signals generated by satellite vehicles (SV&#39;s), civil modes of navigation signals generated by ground and air beacons generated by pseudolites, in the presence of large numbers of malicious tonal and narrowband jammers; malicious signals attempting to spoof those signals; and time and frequency selective multipath induced by natural reflections or malicious repeaters. The system can operate without prior knowledge of, or a search over, the ranging codes for those signals, and even if those signals are received at widely separated power levels; and perform this differentiation based only on the time, frequency, and code diversity between those signals. 
     In one aspect, the multi-subband channelization architecture is extended to spatial and/or polarization diverse multi-feed reception systems, in which each subband employs linear combiners, and linear-algebraic combiner adaptation algorithms, which combine data over spectral and spatial/polarization diversity dimensions to separate GNSS signals and excise jamming interference, and then combine detection features from each subband over the full bandwidth of the GNSS signal. Such systems can exploit spatial and/or polarization diversity between the true navigation signals and the malicious spoofers and repeaters to further blindly detect and differentiate the navigation signals from wideband jammers, and from signals generated by closely aligned repeaters and targeted spoofers that are attempting to match the time, frequency, and code diversity of the navigation signals. As an added advantage, the dimension of each subband combiner can be held constant as the number of spatial/polarization diversity dimensions grows, resulting in a linear growth in complexity as spatial/polarization diversity is increased. 
     Aspects of the disclosure described herein include a system having identifying means comprising at least one antenna, coupled to at least one receiver receptive to energy in a band containing at least one navigation signal with an MOS-DSSS modulation format. The identifying means can include a receiver, for example, dedicated or software-defined radio resources, and a digital signal processor (DSP) (e.g., DSP FPGA, GPU, or the like, with associated storage capability), or combinations of DSP hardware configured for implementing software or firmware installed on that DSP hardware. In these aspects, reception operations can comprise:
         receiving and amplifying at least one signal-in-space (SiS) containing at least a subband of the at least one navigation signal, comprising at least one signal with a MOS-DSSS modulation format comprising a baseband signal with symbol rate f sym , spread by a ranging code that repeats every baseband symbol;   frequency-shifting the at least one SiS to an intermediate or complex baseband frequency, e.g., using a direct-conversion mixer with LO frequency f R  and appropriate filtering operations; and   sampling the frequency-shifted SiS using an analog-to-digital convertor (ADC) with a sampling rate f ADC  that is equal to an integer multiple of the baseband symbol rate f sym , i.e., f ADC =M ADC f sym , where M ADC  is a (typically large) integer.
 
In these aspects, the DSP can perform the following operations:
   capturing at least one sampled and frequency-shifted SiS over a block of received data samples, and recording the approximate reception time and reception frequency of that data block, referred to herein as a “time/frequency stamp”;   channelizing the sampled received data block using at least one symbol-rate synchronous channelizer to provide at least K chn  frequency channels with output rate f sym , each with frequency extent lying within the active bandwidth of the MOS-DSSS navigation signal;   organizing the at least K chn  frequency channels into K sub  subsets of frequency channels, referred to herein as subbands, each comprising M sub  frequency channels;   in multifeed receivers, wherein multiple spatial and/or polarization diverse antennas are coupled to M feed  receivers receptive to energy in a band containing at least one navigation signal with an MOS-DSSS modulation format, channelizing each feed into K chn  frequency channels with output rate f sym , organizing those frequency channels into K sub  subbands each comprising M sub  frequency channels, and combining those subbands across receiver feeds to form a M DoF ×1 vector signal with symbol rate f sym  for each subband, where M DoF =M sub M feed  is referred to herein as the degrees of freedom (DoF&#39;s) of the subband, and where M DoF =M sub  for the single-feed receiver (M feed =1);   collecting the K sub  subbands of M DoF ×1 data vectors over N sym  symbols, and adding a time and frequency stamp those collections of symbols, to form a time/frequency stamped multi-subband output data block;   obtaining third-party baseband symbol data or navigation data from time channels and frequency bands covering the time/frequency stamp, including uncertainty between the receiver LO, clock rate, and/or clock timing phase and universal time coordinates (UTC), using external resources also available to the system; and   processing the time/frequency stamped multi-subband output data block to obtain RPNT analytics, using blind adaptive processing operations responsive to at least one navigation signal received by the system, and using the third-party baseband symbol data or navigation data.
 
In some aspects of the disclosure, the sampled receive data block, or the channelized received data block, along with its time/frequency stamp, might be transported to additional DSPs for processing to obtain the RPNT analytics. These additional DSPs may be onboard the platform containing the reception means, or in a different location entirely, including infrastructure equipment, i.e., “in the Cloud,” if the platform possesses communications means allowing those data blocks to be transported from the platform.
       

     In some aspects, the reception operations are instantiated using purpose-built hardware and DSP hardware. In other aspects, the reception operations are instantiated using general-purpose software-defined radio (GPSDR) and DSP hardware. In either aspect, the receiver clock may be adjusted to remove error from UTC, once a PNT solution is obtained for the receiver. Alternately, the ADC output data may be frequency-sampled and resampled to conform to UTC using algorithms implemented in DSP hardware, once that solution is obtained. 
     In the various aspects, the ADC sampling rate need not have any relation to the chip-rate of the navigation signal; nor need it have any relation to the bandwidth of the navigation signal. In particular, it can have a value that is easily instantiated in low-cost, general-purpose, receiver hardware; and it can have a value that is substantively less than the bandwidth of the navigation signal. 
     In one aspect of the disclosure, the symbol-rate synchronous subband generation is performed by first applying an M ADC -sample fast Fourier transform (FFT) to each S/P output data vector on each of M feed  receiver feeds, and then selecting K FFT  output bins lying within the active bandwidth of the MOS-DSSS navigation signal, to create K chn  frequency channels each with output rate f sym  and dimensionality M feed ×1. These M feed ×1 dimensional frequency channels are then organized into K sub  subbands with M sub  channels per subband, creating K sub  subband signals with output sampling rate f sym  and dimension M DoF ×1, where M DoF =M sub M feed . In another aspect of the disclosure, the subband generation is performed by passing the ADC output signals provided by each of M feed  receiver feeds through a bank of K sub  frequency converters and decimators, each with decimated output sampling rate f dec =M dec f sym , and then passing each decimated signal through a symbol-rate synchronous channelizer to produce a separate M sub ×1 subband signal with sampling rate f sym , on each feed, and then “stacking” those signals to generate a M DoF ×1 vector signal with sampling rate f sym , where M DoF =M sub M feed  However, the channelizer output data may be generated using many different channelizer structures. 
     The DSP software can implement algorithms that exploit baseband symbol sequences or navigation data provided by external or “third-party” resources, for each navigation signal of interest to the system, to efficiently obtain RPNT analytics without the need to demodulate the navigation data. Exemplary RPNT analytics can include:
         estimates of observed received parameters of navigation signals, including signal geo-observables usable to obtain PNT solutions for the reception platform, e.g., observed signal frequency-of-arrival (FOA), time-of-arrival (TOA), and observed local direction-of-arrival (DOA) or global line-of-bearing (LOB) in systems employing multi-feed receivers;   estimates of received signal quality, e.g., received incident power, and of despread/demodulated navigation sequence quality, e.g., despreader output signal-to-interference-and-noise ratio (SINR); and   measures of accuracy of parameter and signal quality estimates.
 
RPNT analytic measurement methods include “partially-blind” methods that do not require knowledge of the ranging code for the navigation signals or spatial/polarization array-manifold data (measurement of cross-feed spatial/polarization signatures as a function of DOA) for the receiver, and “copy-aided” methods that may require one or both of the ranging codes or the array manifold to provide more complete RPNT analytics.
       

     In some aspects, referred to herein as “multi-subband” systems, the blind pilot-exploiting adaptive detection and geo-observable estimation DSP algorithm can be implemented over multiple contiguous or non-contiguous channelizer subbands. Within each subband, the method can detect and separate signals without knowledge of the spreading code; in the presence of arbitrary channel multipath, tonal and narrower-band jammers; and in the presence of near-far interference that is much higher than the spreading gain of conventional PNT systems. The method then combines the subbands to provide a detection statistic and geo-observable estimate that exploits the full bandwidth of the system. The method can also provide a natural means for rejecting narrowband jamming, such as by simply “de-emphasizing” subbands containing excessive interference energy. As part of the FOA estimation step, the method may also compute the differential FOA between subbands, thereby allowing direct estimation of the observed transmitter speed. In aspects where multiple navigation signal transmitter ephemeres are available, the method can provide a full PNT solution based on this information alone. 
     In one aspect, the DSP algorithms can be implemented over multiple portions of the navigation signal band, referred to herein as “navigation signal subbands.” The navigation signal subbands can be obtained using combinations of analog and digital processing. For example, off-center subbands can be obtained by frequency-shifting the SiS with frequency shift f R  that differs substantively from the navigation signal transmit frequency, and sampling the frequency-shifted SiS using an ADC with a sampling rate that is much lower than the bandwidth of the navigation signal. Alternately, off-center subbands can be obtained by downconverting, sampling, and frequency-channelizing the navigation signal over its full bandwidth, and then selecting subsets of frequency channels to for each subband. In both cases, the subbands are processed using the adaptive detection, geo-observable estimation, and (for blind methods) symbol stream estimation DSP algorithms described herein. Any combination of methods spanning this range of methods can be similarly employed. 
     If the number of channelizer DoF&#39;s in each subband is greater than the total number of MOS-DSSS signals (e.g., legitimate navigation signals, spoofers, and substantively time and/or frequency offset multipath images), and tonal interferers in the subband, plus additional processing gain needed to integrate the pilot signal above the non-MOS-DSSS noise-and-interference power level (not limited by the jamming margin of the ranging code), then some aspects of the disclosure can be used to detect and determine geo-observables of all of the MOS-DSSS navigation signals covering that subband without knowledge of the spreading code. Moreover, aspects of the disclosure can be used to perform this detection and geo-observables estimation in the presence of near-far interference that is much higher than the spreading gain of conventional PNT systems, even given the much narrower bandwidth of the single-subband system. In some use scenarios, e.g., reception of navigation signals transmitted from ground or air beacons, where the signal RIP&#39;s are at or near the noise level, little or no additional processing gain is needed to detect the navigation signals, and M DoF  can be much less than the spreading gain of the navigation signals. At the same time, the number of baseband symbols needed to obtain a useful detection solution can be greatly reduced by minimizing M DoF . Consequently, this aspect can sharply reduce the complexity of the DSP algorithms, and the susceptibility of the method to “signature blur” induced by motion of the signal transmitters and receivers over the data collection interval. 
     In some aspects, the navigation signal subbands covered by the channelizer is also dynamically adjusted between collection intervals, e.g., by “stepping” the channelizer through multiple navigation signal subbands, randomly moving the channelizer to different subbands during different collection intervals, or re-tasking the channelizer based on metrics computed during prior collections. This can further reduce susceptibility of the system to narrowband or partial-band jammers attempting to block a portion of the navigation signal band, perhaps a large portion of the band. 
     In one aspect, dynamic adjustment of the channelizer band is used to identify and remove ambiguities in the FOA estimates provided by the system, based on differential Doppler shift between widely-spaced navigation signal subbands. For example, assuming a 1.2 kilometer/second observed SV speed, the GPS L5 short civil signal can experience a 40 Hz difference in Doppler shift over its ˜10 MHz bandwidth (80 Hz if the 20 MHz null-to-null bandwidth is assumed), which is well within the resolution capabilities of the method over 1-to-4 second collect intervals, and is sufficient to estimate the observed SV speed, and thereby resolve the 1 kHz FOA ambiguity induced by the 1 ksps baseband symbol rate of this signal. In aspects where multiple navigation signal transmitter ephemeres (time varying position, velocity, acceleration, and transmission time in universal time coordinates) are available, e.g., determined from demodulated navigation data or external resources, methods can provide a full PNT solution based on this information alone. 
     The DSP detection and geo-observable estimation algorithm can extend without qualitative change to scenarios in which the receiver possesses M feed ≥1 spatial/polarization diverse receive feeds, except that the channelizer DoF&#39;s are multiplied by the number of receive feeds, such that M DoF =M chn ·M feed , where M chn  is the number of channels per feed. However, the resultant system can additionally exploit spatial/polarization diversity of MOS-DSSS signals, including target spoofers that closely emulate the time, frequency, and code diversity of the legitimate navigation signals, and can further excise up to M feed −1 wideband jammers, on the basis of said spatial/polarization diversity, even if those jammers are much stronger than the legitimate navigation signals. As in single-feed aspects, the total DoF&#39;s must additionally be greater than at least the number of MOS-DSSS signals (e.g., legitimate navigation signals and spoofers) and tonal signals impinging on the receiver, over the bandwidth used by the system. 
     In one aspect of the disclosure, M chn ≈M DoF /M feed , such that M DoF  is held constant as the number of feeds is increased. This can further reduce susceptibility of the system to signature blur and provide enhanced rejection of jammers and targeted spoofers attempting to match the time, frequency, and code diversity of legitimate navigation signals, even if those jammers or spoofers are much stronger than the legitimate navigation signals. 
     In another aspect of the disclosure, a multi-subband channelizer with K sub  subbands is implemented on each receiver feed, and M chn  and K sub  are simultaneously set to M chn ≈M DoF /M feed  and K sub ≈M DoF ·M feed , such that both M DoF  and the total bandwidth used by the system are held constant as the number of feeds is grown. Assuming the constraints on M DoF  and M feed  given above, the resultant system can achieve the full processing gain of the navigation signal modulation format, and the robustness to near-far interference exhibited by the method in any reception scenario, without unduly affecting complexity or convergence time of the system. 
     In some aspects of the disclosure, copy-aided post-processing geo-observable estimation methods, along with ranging code information determined from demodulated navigation data or through external resources, are used to resolve the FOA and TOA ambiguities. In systems employing multi-feed receivers, further aspects can use the directional array manifold of the receiver to determine the direction of arrival (DOA) of the received signals; the receiver orientation; and true lines-of-bearing (LOB&#39;s) from the receiver to the signal emitters. In both aspects, the fine geo-observables can be used to geolocate spoofers identified by the system, using known location of the receiver and geo-observables of those spoofers. 
     In another aspect of the disclosure, the system employs additional network communication methods to allow intercommunication of parameters between multiple receivers. These receivers can employ combinations of single-feed or multi-feed receivers, and can have different ADC rates and channelization methods, as long as each channelizer provides a vector sequence at a 1 ksps output data rate. Among other advantages, this aspect provides additional means for detecting and defeating target spoofers, by using additional network diversity to prevent targeted spoofers from overwhelming any single receiver in the environment. 
     Any of these aspects can demodulate civil navigation signals in the presence of arbitrary multipath because they do not explicitly exploit the PRN code in the demodulation process. Any of these aspects can provide for determining the yaw-pitch-roll orientation of the receive platform, such as due to relationship between Doppler shift and direction vectors to the GNSS emitters, once the GNSS emitter and receiver positions are determined. 
     In networked receiver configurations, intercommunication can allow each receiver to continuously draw from and update a “network database” containing estimated TOA&#39;s, FOA&#39;s, and (for multi-feed receivers) DOA&#39;s for navigation signals detected by any receiver in the network, and containing ephemeris information for those receivers. This database can be used to further identify spoofers, including targeted spoofers, based on differential FOA (DFOA) and TOA (DTOA) measured at those receivers, which removes any deliberate timing or frequency shift used by the spoofer to emulate a true navigation signal at a specific victim receiver. If the victim receiver does not employ a multi-feed receiver, and may therefore find itself being successfully spoofed by a targeted attack, this network will provide an unambiguous alert to such an attack, as the spoofed signals will not possess the true GNSS FOA or TOA at the other receivers in the network, and will therefore be removed using the blind adaptive despreading methods employed in aspects of the disclosure. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
       Block diagrams and flow diagrams depicted herein can represent computer software instructions or groups of instructions. One or more of the processing blocks or steps may be implemented with modules (e.g., procedures, functions, and so on) that perform the functions described herein. The software codes may be stored in non-transitory computer-readable memory. Alternatively, processing blocks or steps described herein may represent steps performed by hardware, including one or more functionally equivalent circuits, such as a digital signal processor, an application specific integrated circuit, a programmable logic device, a field programmable gate array, a general-purpose processor programmed to perform one or more of the disclosed steps, or other electronic units designed to perform the functions disclosed herein. Disclosed aspects can employ multiple processors communicatively coupled together, and can employ distributed computing, such as Cloud computing, cluster computing, distributed computing, etc. 
         FIG.  1    depicts a system in which aspects of the disclosure can be implemented. 
         FIG.  2    illustrates operating parameters and geo-observables that can be employed in aspects disclosed herein. 
         FIG.  3    is a block diagram that depicts generation of civil navigation signals that aspects of the disclosure can be configured to exploit. 
         FIG.  4    is a Table of civil GNSS modulation formats that aspects of the disclosure can be configured to exploit. 
         FIG.  5 A ,  FIG.  5 B ,  FIG.  6 A , and  FIG.  6 B  are block diagrams of receiver front-ends that can be configured according to aspects of the disclosure. 
         FIG.  7    is a block diagram that depicts a single-feed channelizer in accordance with an aspect of the disclosure. 
         FIG.  8    is a block diagram that depicts a single-feed FFT-based channelizer in accordance with an aspect of the disclosure. 
         FIG.  9    is a block diagram of a single-feed FFT-based symbol-rate synchronous multi-subband channelizer according to an aspect of the disclosure. 
         FIG.  10    is a block diagram of a single-feed polyphase filter based symbol-rate synchronous multi-subband channelizer according to an aspect of the disclosure. 
         FIG.  11    is a block diagram of single-feed symbol-rate synchronous multi-subband channelizer according to another aspect of the disclosure. 
         FIG.  12 A  and  FIG.  12 B  depict exemplary channelizer metrics for a rectangularly windowed polyphase filter based symbol-rate synchronous channelizer. 
         FIG.  13    is a block diagram that depicts a multifeed receiver and channelizer in accordance with an aspect of the disclosure. 
         FIG.  14    is a block diagram of a multifeed FFT-based symbol-rate synchronous multi-subband channelizer according to an aspect of the disclosure. 
         FIG.  15    is a block diagram of a multifeed polyphase filter based symbol-rate synchronous multi-subband channelizer according to an aspect of the disclosure. 
         FIG.  16    and  FIG.  17    illustrate despreader parameters that can be employed in aspects disclosed herein. 
         FIG.  18    is a flow diagram of a multi-subband signal detection, geo-observable/quality estimation, and symbol estimation algorithm that can be employed in aspects disclosed herein. 
         FIG.  19    is a flow diagram of a multi-subband signal detection and geo-observable/quality estimation, algorithm that can be employed in aspects disclosed herein. 
         FIG.  20    is a flow diagram of a multi-subband signal detection and geo-observable/quality estimation, algorithm that can be employed in aspects disclosed herein. 
         FIG.  21    is a flow diagram of a local-maxima search procedure according to an aspect of the disclosure. 
         FIGS.  22 A,  22 B, and  22 C  and  FIGS.  23 A,  23 B, and  23 C   FIG.  22    and  FIG.  23    are plots that depict DFOA search sensitivity. 
     
    
    
     DETAILED DESCRIPTION 
     Various aspects of the disclosure are described below. It should be apparent that the teachings herein may be instantiated in a wide variety of forms and that any specific structure, function, or both being disclosed herein are merely representative. Based on the teachings herein one skilled in the art should appreciate that an aspect disclosed herein may be implemented independently of any other aspects and that two or more of these aspects may be combined in various ways. For example, an apparatus may be implemented or a method may be practiced using any number of the aspects set forth herein. In addition, such an apparatus may be implemented or such a method may be practiced using other structure, functionality, or structure and functionality in addition to or other than one or more of the aspects set forth herein. 
     In the following description, for the purposes of explanation, numerous specific details are set forth in order to provide a thorough understanding of the invention. It should be understood, however, that the particular aspects shown and described herein are not intended to limit the invention to any particular form, but rather, the invention is to cover all modifications, equivalents, and alternatives falling within the scope of the invention as defined by the claims. 
       FIG.  1    depicts a model of the presumed normal condition wherein a receiver ( 1 ) is receiving GPS signals from an SV ( 3 ). 
       FIG.  2    Illustrates 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 ( 6 ) 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 (CCP) (e.g., after calibration and removal of transmitter and receiver phase offsets), which are measurable and differentiable at the receiver ( 1 ). 
     The TOA, RIP, FOA, LOB, DOA, and CCP 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 positioning, 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. 
       FIG.  3    depicts operations used to generate the civil navigation signals that are the focus of the invention. It should be noted that the processing steps shown in  FIG.  3    is one conceptual means for generating civil GNSS signals, and ignore other SiS&#39;s that may be transmitted simultaneously with civil GNSS signals, for example, the GPS P, P(Y), M, and L1C SiS&#39;s. In practice, these other SiS&#39;s are largely outside the band occupied by the civil GNSS signals, and moreover are received well below the noise floor. Therefore, their effect on aspects 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. 
     The processing steps shown in  FIG.  3    generate a particular class of spread spectrum signal, referred to herein as a modulation-on-symbol direct-sequence spread spectrum (MOS-DSSS) signal, in which a data symbol sequence d T (n sym ;l) with symbol rate f sym  is spread by a M chp ×1 dimensioned ranging code 
                 c   T     (   ℓ   )     =       [       c   T     (       m   chp     ;   ℓ     )     ]         m   chp     =   0         M   chp     -   1             
that is repeated over every data symbol. The spread signal can be represented as the multiplication of data symbol sequence d T (n sym ;l) and code vector c T (l) ( 115 ), resulting in M chp ×1 c T  (l)d T (n sym ;l). The vector signal is then passed through a M chp :1 parallel-to-serial (P/S) converter ( 116 ), resulting in a scalar spread signal with rate M chp f sym , which can be expressed as c T (n chp  mod M chp ;l)d T (└n chp /M chp ┘;l) where (⋅)mod M and └⋅┘ denote the modulo-M operation and integer truncation or “floor” operations, respectively. As shall be described below, this property induces inherent massive spectral redundancy to the civil GNSS format that is exploited by disclosed aspects.
 
     The chip-rate spread signal is then passed through digital-to-analog conversion (DAC) ( 118 ), and transmission/upconversion stages ( 120 ), controlled by an internal clock (clk) ( 117 ) and local-oscillator (LO) ( 119 ) stages, and transmitted from an antenna ( 121 ) to generate a signal-in-space (SiS) given by √{square root over (2P T (l))}Re{s T (t−τ T (l);l)e j(2πd     T     t+φ     T     (l)) } for transmitter l ( 122 ), where f T (l) is the transmit frequency, typically common to all SV&#39;s (f T (l)≡f T ), P T (l), φ T (l) and τ T (l) are the transmit power, phase, and electronics delay for SiS l, respectively, and where s T (t;l) is the complex baseband representation of the transmitted signal. 
     Based on these operations, the complex-baseband signal can be modeled by 
                         s   T     (     t   ;   ℓ     )     =       ∑           n   sym               d   T     (       n   sym     ;   ℓ     )     ⁢       h   T     (       t   -       T   sym     ⁢     n   sym         ;   ℓ     )           ,           Eqn   ⁢           (   1   )                 
where d T (n sym ;l) is the encoded symbol stream, T sym =1/f sym , is the symbol period, and h T (t;l) is a spread symbol given by
 
                               h   T     ⁢     (     t   ;   ℓ     )       =         ∑       m   chp     =   0           M   chp     -   1           c   T     ⁢     (       m   chp     ;   ℓ     )     ⁢     g   T     ⁢     (       t   -       T   chp     ⁢     m   chp         ;   ℓ     )           ,               T   chp     =       T   sym       M   chp         ,                 Eqn   ⁢           (   2   )                 
and where g T  (t;l) is the chip symbol modulated by each chip, referred to here as the chip shaping. The spread symbol can also be expressed by it&#39;s continuous Fourier transform H T (f;l)=∫h T (f;l)e −j2πft  df,
 
 H   T ( f;l )= C   T ( e   j2πfT     chp     ;l ) G   T ( f;l ),  Eqn (3)
 
where
 
                 C   T     (       e     j   ⁢   2   ⁢   π   ⁢   f       ;   ℓ     )     =       ∑     m   chp             c   T     (     m   chp     )     ⁢     e       -   j     ⁢   2   ⁢   π   ⁢     fm   chp                   
is the discrete Fourier transform of C T (m chp ;l), and where G T (f;l) is the continuous Fourier transform of g T (t;l). The spread signal can therefore be modeled as a pulse-amplitude modulated (PAM) waveform in which baseband symbol sequence modulates a symbol waveform that is itself direct-sequence spread by the ranging code, leading to the MOS-DSSS nomenclature employed here.
 
     In typical navigation systems, the chip shaping is identical for each transmitted signal (g T (t;l)≡g T (t)); however, in practice, the shaping can also vary between transmitters. For typical navigation systems, the chip shaping is also nominally rectangular, such that G T (f;l)≡T chp Sa(πfT chp ), where Sa(x)=sin(x)/x; however, different chip shapings can be employed without changing the basic structure of the MOS-DSSS modulation format. 
     In typical navigation systems, the spreading code C T (l) employed at transmitter l is unique to that transmitter, and is in fact provisioned by the navigation system. In contrast, the transmit power and center frequency are typically (but not necessarily) identical for each transmitter, while the transmitter carrier phase, delay, and other impairments will vary between the transmitters. However, any of these typical assumptions can be varied in practice, or in different instantiations of a navigation system. 
     As shown in  FIG.  3   , the full symbol sequence d T (n sym ;l) is generated from an underlying navigation sequence b T (n NAV ;l)| NAV  with navigation rate f NAV , which contains information needed to generate a full positioning, navigation, and timing (PNT) solution at navigation receiver. The navigation sequence typically provides timing and ephemeris (time-stamped position, velocity, and acceleration) of the navigation transmitter and larger network, known impairments at the transmitter, and additional information needed to determine quality and timing of the transmitted signal. The navigation sequence also typically provides information about the specific spreading code(s) employed at the transmitter. This navigation sequence is then passed through processing operations to generate an encoded signal sequence b T (n sym ;l)| sym , with symbol rate f sym  and for some modulation formats combined with a known pilot sequence p T (n sym ;l), with symbol rate f sym  ( 112   b ) to generate the full symbol sequence signal d T (n sym ;l). All of these additional processing steps add additional known structure to the transmitted MOS-DSSS signal, which can be exploited in subsequent processing stages. Additionally, the navigation sequence, the information contained within that sequence (e.g., the transmitter ephemeres) can be obtained after transmission from third-party sources, providing additional exploitable structure in applications some aspects of the invention. 
       FIG.  4    is a tabular description of GNSS networks employing civil GNSS signals that can be modeled as 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.  3   . 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. 
     As the last column of  FIG.  4    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) or “legacy” signal listed in the first row of  FIG.  4    is constructed from a 50 bit/second (bps) BPSK (real) navigation sequence b T (n NAV ;l)| NAV , repeated 20 times per MOS-DSSS symbol to generate 1 ksps symbol sequence d T (n sym ;l)=b T (n sym ;l)| sym =b T (└n sym /M NAV ┘;l)| NAV  where M NAV =20. Moreover, b T (n NAV ;l)| NAV  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. 
     Similarly, the GPS L5, QZSS, Galileo E5B, E6B, and Beidou 1B transmitters all add a periodic pilot to the quadrature rail of their navigation signals, and the NAVIC RS BOC transmitter adds a periodic pilot to the in-phase rail of its navigation signal. Means for exploiting periodic pilots are a focus of some aspects of this invention. 
       FIG.  5 A  and  FIG.  5 B  each depict exemplary direct-conversion front-end reception processing for a single-feed receiver receiving GNSS signals.  FIG.  5 A  shows a direct-conversion receiver front-end that is dedicated to a specific GNSS band, and  FIG.  5 B  shows 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, some of the disclosed aspects. 
     The direct-conversion receivers shown in  FIG.  5 A  and  FIG.  5 B  each comprise the following:
         An antenna ( 128 ) and low-noise amplifier (LNA) ( 131 ), which receives an SiS containing GNSS emissions within the FoV of the receiver.   A local oscillator (LO) ( 133 ), which generates a complex sinusoid with frequency f R  and phase φ R , comprising an in-phase (I) real sinusoid, and a quadrature (Q) sinusoid shifted by 90° in phase from the in-phase sinusoid;   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 R      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 R (t), centered at 0 MHz if f R =f T .   A dual analog-to-digital converter (ADC) ( 136 ), which samples each rail of x R (t) at rate f ADC =M ADC f sym , where f sym  is the civil GNSS signal symbol rate shown in  FIG.  3    and M ADC  is a positive integer, thereby generating complex dual-ADC output signal x ADC (n ADC ) with time index n ADC  ( 137 ).
 
The receiver structure shown in  FIG.  5 A  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.
       

     In contrast, the receiver structure shown in  FIG.  5 B  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 R  is significantly offset from the GNSS band center f 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 can be an advantageous feature of at least some of the disclosed aspects. 
       FIG.  6 A  and  FIG.  6 B  each depict exemplary front-end reception processing for a single-feed superheterodyne receiver receiving GNSS signals; with  FIG.  6 A  showing a superheterodyne receiver front-end that is dedicated to a specific GNSS band, and  FIG.  6 B  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 some aspects disclosed herein. 
     The heterodyne receivers shown in  FIG.  6 A  and  FIG.  6 B  shift the SiS received by an antenna ( 128 ) and amplified in an LNA ( 130 ) to IF frequency f IF , typically using a two-stage mixer  134 ,  141 , yielding analog scalar real-IF signal x R (t) centered on IF frequency f IF  if f R =f T , with receive phase-shift φ R =φ R   (1) +φ R   (2) . The analog signal is then passed to an ADC ( 143 ) that samples x R (t) at rate M ADC f sym , where f sym  is the MOS-DSSS signal symbol rate shown in  FIG.  6 A  and  FIG.  6 B , and M ADC  is a positive integer, yielding real ADC output signal X ADC (n ADC ) with time index n ADC  ( 137 ). The receiver structure shown in  FIG.  6 A  can be optimized for a specific GNSS band, while the receiver structure shown in  FIG.  6 B  can be flexibly tasked to different navigation signal bands, or to segments of such bands, under software control. 
     In all four cases, the receiver can induce additional delays, frequency offsets, and sampling rate offsets due to internal electronics in the receiver systems, and due to clock rate and timing offset from universal time coordinates (UTC). These are subsumed into the overall channel response, and measured and removed as part of the geo-location solution. 
     As will be discussed below, the ADC output signal provided by either the direct-conversion or the superheterodyne receivers can be used in disclosed aspects without any change to the design, except control parameters used in channelization operations immediately after that ADC. The disclosed aspects can also be practiced with many other receiver front-ends known to the art. In many cases, the disclosed aspects are blind to substantive differences between those receiver front-ends. 
     Assuming the reception scenario shown in  FIG.  2   , and the direct-conversion receiver structure shown in  FIG.  5   , the downconverted and filtered signal shown in x R (t) is given by 
                         x   R     (   t   )     =         i   R     (   t   )     +         ∑     ℓ   =   1       L         s   R     (     t   ;   ℓ     )           ,           Eqn   ⁢           (   4   )                 
where i R (t) is the noise and interference received in the receiver passband and s R (t;l) is the complex-baseband representation of the signal received from transmitter l. Over time intervals on the order of 1-10 seconds from reception time t 0 , s R (t;l) can be modeled by
 
                         s   R     (         t   0     +   t     ;   ℓ     )     ≈         g   TR     (   ℓ   )     ⁢     e     j   ⁢   2   ⁢     π   ⁡   (           α   TR     (   ℓ   )     ⁢   t     +       1   2     ⁢   α   ⁢           (   1   )             TR         ⁢     (   ℓ   )     ⁢     t   2         )         ⁢       s   T     (           (     1   +     τ   ⁢           (   1   )             TR         ⁢     (   ℓ   )         )     ⁢   t     -       τ   TR     (   ℓ   )       ;   ℓ     )         ,           Eqn   ⁢           (   5   )                 
at the input to the dual-ADC&#39;s shown in  FIG.  5   , where
 
                         g   TR     (   ℓ   )       =   △             P   R     ⁢     (       t   0     ;   ℓ     )         ⁢     exp   ⁡   (     j   ⁢       φ   TR     (   ℓ   )       )         ,           end   ‐   to   ‐   end   ⁢         complex   ⁢         channel   ⁢         gain             Eqn   ⁢           (   6   )                                   φ   TR     (   ℓ   )       =   △           φ   T     (   ℓ   )     -     φ   R     -     2   ⁢   π   ⁢       f   T     (   ℓ   )     ⁢       τ   TR     (       t   0     ;   ℓ     )           ,           end   ‐   to     ⁢   ‐   ⁢     end   ⁢         channel   ⁢         phase             Eqn   ⁢           (   7   )                                         τ   TR     (   ℓ   )       =   △           τ   T     (   ℓ   )     +     τ   R     +       τ   TR       (   0   )             (       t   0     ;   ℓ     )         ,           end   ‐   to   ‐   end   ⁢         observed   ⁢                ⁢   T   ⁢   O   ⁢   A                 Eqn   ⁢           (   8   )                                         α   TR     (   ℓ   )       =   △           f   T     (   ℓ   )     +     f   R     -         f   T     (   ℓ   )     ⁢       τ   TR     (   1   )       (       t   0     ;   ℓ     )           ,           end   ‐   to   ‐   end   ⁢         observed   ⁢         F   ⁢   O   ⁢   A                 Eqn   ⁢           (   9   )                                         τ   TR     (   1   )       (   ℓ   )       =   △         τ   TR     (   1   )       (       t   0     ;   ℓ     )       ,           observed   ⁢         differential   ⁢         T   ⁢   O   ⁢     A   ⁡   (     D   ⁢   T   ⁢   O   ⁢   A     )                   Eqn   ⁢           (   10   )                                         α   TR     (   1   )       (   ℓ   )       =   △         -       f   T     (   ℓ   )       ⁢       τ   TR     (   2   )       (       t   0     ;   ℓ     )         ,           observed   ⁢         differential   ⁢         F   ⁢   O   ⁢       A   ⁡   (     D   ⁢   F   ⁢   O   ⁢   A     )     .                   Eqn   ⁢           (   11   )                 
and where φ T (l), f T (l), and τ T (l) are the carrier phase, frequency, and timing offset induced at transmitter l; φ R , f R , and τ R  are the carrier phase, frequency, and timing offset induced at the receiver; P R (t 0 ;l) is the received incident power of SiS l at the receiver; and τ TR   (q) (t 0 ;l)=dτ TR (t;l)/dt q  is the q th  derivative of the observed TOA between transmitter l and the receiver at time t. Given observed position, velocity, and acceleration p TR (t;l)=p T (t;l)−p R (t;l), v TR (t;l)=v T (t;l)−v R (t;l), and a TR (t;l)=a T (t;l)−a R (t;l), where {p T (t;l),v T (t;l),a T (t;l)} l=1   L  and {p R (t),v R (t),a R (t)} are the transmitter and receiver ephemeres, respectively, then τ TR   (q) (t;l) is given by
 
                         τ   TR     (   0   )       (     t   ;   ℓ     )     =       1   c     ⁢              p   TR     (     t   ;   ℓ     )          2         ,           Eqn   ⁢           (   12   )                                   τ   TR     (   1   )       (     t   ;   ℓ     )     =       1   c     ⁢       u   TR   T     (     t   ;   ℓ     )     ⁢       v   TR     (     t   ;   ℓ     )         ,           Eqn   ⁢           (   13   )                                   τ   TR     (   2   )       (     t   ;   ℓ     )     =       1   c     ⁢     (           u   TR   T     (     t   ;   ℓ     )     ⁢       a   TR     (     t   ;   ℓ     )       +                P   ⊥     (       p   TR     (     t   ;   ℓ     )     )     ⁢       v   TR     (     t   ;   ℓ     )            2       )         ,           Eqn   ⁢           (   14   )                 
where u TR (t;l) is the line-of-bearing (LOB) direction vector from the receiver to transmitter l,
 
                         u   TR     (     t   ;   ℓ     )     =         p   TR     (     t   ;   ℓ     )                p   TR     (     t   ;   ℓ     )          2         ,           Eqn   ⁢           (   15   )                 
and P ⊥ (p TR (t;l)) is the projection matrix orthogonal to p TR (t;l), given by
 
                         P   ⊥     (       p   TR     (     t   ;   ℓ     )     )     =       I   3     -           p   TR     (     t   ;   ℓ     )     ⁢       p   TR   T     (     t   ;   ℓ     )                  p   TR     (     t   ;   ℓ     )          2   2           ,           Eqn   ⁢           (   16   )                 
and where ∥⋅∥ 2  and I 3  denote the L2 Euclidean vector norm and 3×3 identity matrix, respectively. As shown in the &#39;417 NPA, incorporated herein by reference, the RIP and LOB adheres closely to a zero-order model over time intervals on the order of 1-to-10 seconds, for transmitters in medium-Earth orbits and receivers with dynamics commensurate with high-velocity airborne vehicles.
 
     The local direction-of-arrival (DOA) of the signal received from transmitter l can then be represented by local DOA direction vector u R (t;l)=Ψ R (t)u 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. As shown in the &#39;417 NPA, incorporated herein by reference, the DOA also adheres closely to a first-order model over time intervals on the order of 1-to-10 seconds, for transmitters in medium-Earth orbits and receivers with dynamics commensurate with high-velocity airborne vehicles. Moreover, the variation in local DOA is nearly identical for all of the received navigation signals under this scenario, as it primarily due to motion/orientation of the receiver, which equally affects all of the signal LOB&#39;s. 
       FIG.  7    shows a single-feed receiver structure used in one aspect of the disclosure, describing the symbol-rate synchronous vector channelizer ( 150 ) integral to the disclosure. 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.  5 A  and  FIG.  5 B , comprising reception, frequency downconversion by receive frequency f R  to complex baseband representation, and sampling using a dual ADC ( 136 ) with sampling rate M ADC f sym , where M ADC  is a positive integer. The sampled signal is then passed to a symbol-rate synchronous channelizer ( 150 ) comprising the following conceptual operations:
         A 1:M ADC  serial-to-parallel conversion operation that transforms the rate M ADC f sym  ADC output signal to an M ADC ×1 vector signal x sym (n sym ) with rate f sym , given by       

     
       
         
           
             
               
                 
                   
                     
                       x 
                       sym 
                     
                     ( 
                     
                       n 
                       sym 
                     
                     ) 
                   
                   = 
                   
                     
                       ( 
                       
                         
                           
                             
                               
                                 x 
                                 R 
                               
                               ( 
                               
                                 
                                   T 
                                   sym 
                                 
                                 ⁢ 
                                 
                                   n 
                                   sym 
                                 
                               
                               ) 
                             
                           
                         
                         
                           
                             ⋮ 
                           
                         
                         
                           
                             
                               
                                 x 
                                 R 
                               
                               ( 
                               
                                 
                                   T 
                                   sym 
                                 
                                 ( 
                                 
                                   
                                     n 
                                     
                                       sym 
                                         
                                     
                                   
                                   + 
                                   
                                     
                                       m 
                                       ADC 
                                       
                                         - 
                                         1 
                                       
                                     
                                     
                                       M 
                                       ADC 
                                     
                                   
                                 
                                 ) 
                               
                               ) 
                             
                           
                         
                       
                       ) 
                     
                     . 
                   
                 
               
               
                 
                   Eqn 
                   ⁢ 
                       
                   
                     ( 
                     17 
                     ) 
                   
                 
               
             
           
         
       
         
         
           
             An M ADC :K chn  channelization operation ( 155 ), conceptually performed by operating on x sym (n sym ) using K chn ×M ADC  channelization matrix operator T chn  ( 153 ) to form K chn ×1 channelizer output signal x chn (n sym )=T chn ∘x sym (n sym ) ( 159 ), where K chn  is the symbol-rate synchronous output channels used in subsequent DSP operations and “∘” denotes a general filtering operation. In some aspects of the invention, T chn  ( 153 ) is a memoryless matrix and the operation is a simple matrix multiplication. In other aspect of the invention, T chn  ( 153 ) is a more linear-time-invariant (LTI) impulse response, and the operation is a matrix convolution.
 
The data rate of x chn (n sym ) is nominally equal to the symbol rate of the navigation 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.  4   , and the dimensionality of x chn (n sym ) is K chn ×1, where K chn  Is the number of frequency channels used by subsequent DSP processing stages. 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. 
       FIG.  8    depicts a specific instantiation of the symbol-rate synchronous vector channelizer ( 150 ), using Fast-Fourier Transform (FFT) operations. The channelizer first applies a (nominally zero-padded) discrete Fourier transform (DFT) with M ADC  input samples, K FFT ≥M ADC  output bins ( 152 ), and a (nominally rectangular) window 
                     {       w   chn     (     m   ADC     )     }         m   ADC     =   0         M   ADC     -   1             (   156   )               
to each S/P output vector, where K FFT  can be conveniently chosen to minimize computational complexity and/or cost of the DFT using an FFT algorithm. The channelizer then selects K chn  output bins
 
               {       k   bin     (     (     k   chn     )       }         k   chn     =   0         K   chn     -   1           
for use in subsequent adaptive processing algorithms ( 154 ). The channelization operation can be expressed as a multiplication of x sym (n sym ) by K chn ×M ADC  matrix
 
     
       
         
           
             
               
                 
                   
                     T 
                     chn 
                   
                   = 
                   
                     
                       [ 
                       
                         
                           exp 
                           ⁡ 
                           ( 
                           
                             
                               - 
                               j 
                             
                             ⁢ 
                             2 
                             ⁢ 
                             π 
                             ⁢ 
                             
                               
                                 k 
                                 bin 
                               
                               ( 
                               
                                 k 
                                 chn 
                               
                               ) 
                             
                             ⁢ 
                             
                               
                                 m 
                                 ADC 
                               
                               
                                 K 
                                 FFT 
                               
                             
                           
                           ) 
                         
                         ⁢ 
                         
                           
                             w 
                             chn 
                           
                           ( 
                           
                             m 
                             ADC 
                           
                           ) 
                         
                       
                       ] 
                     
                     . 
                   
                 
               
               
                 
                   ( 
                   153 
                   ) 
                 
               
             
           
         
       
     
     In one aspect of the invention, the output channelizer bins are selected contiguously over the passband of the receiver, e.g., using formula k bin (k chn )=(k init +k chn )mod K FFT , where k init  is the first FFT bin output from the channelizer (FFT bin associated with channelizer bin 0). As will be discussed herein, an advantage of this approach is minimization of channel signature blur due to movement of the transmitter and/or receiver platforms. In another aspect of the invention, the output channelizer bins are selected sparsely over the passband of the receiver, e.g., using formula k bin (k chn )=(k init +K sep k chn )mod K FFT , where K sep  is an integer separation factor between channelizer bins. An advantage of this approach is ability to avoid known receiver impairments such as LO leakage (e.g., by setting k init  to a non-multiple of K sep ), or to reduce complexity of the channelization operation, e.g., using sparse FFT methods. 
     In some aspects of the invention, 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. In other aspects of the invention, the output channelizer bins are adaptively selected, e.g., to avoid dynamic narrowband co-channel interferers (NBCCI). 
     In, the &#39;417 NPA, which has been incorporated herein by reference, means for processing the vector channelizer output signal, in that case with K chn =M DoF . However, exploitation of channelization dimensionalities covering the full active bandwidth of the GPS L1 legacy signal, much less the much wider bandwidth GPS L5 and other civil GNSS signals, would substantively increase both the complexity of the linear-algebraic methods disclosed in the &#39;417 NPA, and the reception time needed to implement those methods. For this reason, the &#39;417 NPA focused on partial-band and subband methods that only exploited a part of that active bandwidth. Methods for overcoming this limitation without unduly affecting complexity or reception time is a primary focus of this invention. 
       FIG.  9    depicts an instantiation of a single-feed FFT-based symbol-rate synchronous multi-subband channelizer, intended to overcome these limitations. This instantiation is identical to the FFT-based symbol-rate synchronous channelizer shown in  FIG.  8   , except that the K chn -channel selection operation ( 154 ) is replaced by a subband selection operation ( 157 ) that selects K sub  subbands from the K FFT  FFT output bins provided by the FFT operation ( 152 ). In the aspect of the invention shown in this Figure, the number of bins in each subband is equal to the same number, M DoF , referred to herein as the number of degrees of freedom of the subband. This has the benefit of regularizing linear-algebraic DSP operations performed on each subband in subsequent processing steps. However, other aspects of the invention can vary the number of bins, and therefore processing degrees of freedom, e.g., based on numbers of interferers present in each subband or other factors. The FFT output bins selected for the subbands are then organized into K sub  M DoF ×1 vector signals given by 
                         x   DoF     (       n   sym     ;     k   sub       )     =     (             x   FFT     (         k   bin     (     0   ;     k   sub       )     ,     n   sym       )             ⋮               x   FFT     ⁢     (         k   bin     (         M   DoF     -   1     ;     k   sub       )     ,     n   sym       )             )       ⁢   
         k   sub     =   0     ,   …         ,       K   sub     -   1     ,             Eqn   ⁢           (   18   )                 
each with output rate f sym , where
 
               {       x   FFT     (       k   bin     ,     n   sym       )     }         k   bin     =   0         K   FFT     -   1           
are the bins output from the FFT ( 152 ).
 
     In one aspect of the invention, the FFT bins selected for each subband are contiguous, e.g., set using formula k bin (m DoF ;k sub )=(k init (k sub )+m DoF )mod K FFT , where k init (k sub ) denotes the initial FFT bin selected for subband k sub . As will be discussed herein, an advantage of this approach is minimization of channel signature blur due to movement of the transmitter and/or receiver platforms. In another aspect of the invention, the FFT bins selected for each subband are more sparsely distributed, e.g., using formula k bin (m DoF ;k sub )=(k init (k sub )+K sep m DoF )mod K FFT , where K sep  is an integer separation between bins within the subband. An advantage of this approach is ability to avoid known receiver impairments such as LO leakage, or to reduce complexity of the channelization operation, e.g., using sparse FFT methods. The bin spacing K sep  can be equal for each subband, or be different between subbands. 
     In one aspect of the invention, the subband are themselves contiguous, e.g., such that k init (k sub +1)=(k init (k sub )+M DoF )mod K FFT  for subbands with contiguous bins within each subband, and k init (k sub +1)=(k init (k sub )+K sep M DoF )mod K FFT  for subbands with sparsely separated bins within each subband. In other aspects of the invention, subbands may be widely separated. In yet other aspects of the invention, subbands may be interleaved, e.g., by setting K sep =K sub  and k init (k sub )=k init (0)+k sub , such that k bin (m DoF ;k sub )=(k init (0)+(K sub m DoF +k sub ))mod K FFT . 
     In some aspects of the invention, the output subband bins are preselected, e.g., to avoid known receiver impairments such as LO leakage, or to reduce complexity of the channelization operation. In other aspects of the invention, the output subband bins are adaptively selected, e.g., to avoid dynamic narrowband co-channel interferers (NBCCI). 
       FIG.  10    is a block diagram of an alternate single-feed polyphase filter based symbol-rate synchronous multi-subband channelizer, according to an aspect of the disclosure. This block diagram also depicts a simple means for separating the channelizer output signal into contiguous subbands. The ADC output data x R (T ADC n ADC ) is passed through a 1:M ADC  serial-to-parallel converter ( 151 ) and a Q ADC M ADC -tap polyphase filter ( 160 ) with (typically real, potentially complex) channelizer weights 
                       {       w   chn     (     m   ADC     )     }         m   ADC     =   0         Q   ADC     ⁢     M   ADC         ,         158             
given by
 
                         x   chn     (       k   chn     ,     n   sym       )     =       ∑       m   ADC     =   0           Q   ADC     ⁢     M   ADC       -   1             w   chn   *     (     m   ADC     )     ⁢       x   R     (       T   ADC     (         M   ADC     ⁢     n   sym       +     m   ADC       )     )     ⁢     e       -   j     ⁢   2   ⁢   π   ⁢       f   chn     (     k   chn     )     ⁢     T   ADC     ⁢     m   ADC               ,           Eqn   ⁢           (   19   )                 
such that x chn (n sym ) is given by
 
                                              x   chn     (     n   sym     )     =         ∑       Q   ADC     -   1           q   ADC     =   0             T   chn     (     q   sym     )     ⁢       x   sym     (       n   sym     +     q   sym       )           ⁢             T   chn     (     q   sym     )     =           [         w   chn   *     (         M   ADC     ⁢     q   sym       +     m   ADC       )     ⁢     e       -   j     ⁢   2   ⁢   π   ⁢       f   chn     (     k   chn     )     ⁢       T   ADC     (         M   ADC     ⁢     q   sym       +     m   ADC       )           ]     .                 Eqn   ⁢           (   20   )                 
Relating to  FIG.  7   , this channelizer employs a transformation x chn (n sym )=T chn ∘x sym (n sym ), where T chn  is a K chn ×M ADC  polyphase filter with response length Q ADC . In the aspect shown in  FIG.  10   , the channel frequencies
 
               {       f   chn     (     k   chn     )     }         k   chn     =   0         K   chn     -   1           
are furthermore assumed to satisfy
 
                       f   chn     (     k   chn     )     =       (       k   chn     -         K   chn     -   1     2       )     ⁢     f   sym               Eqn   ⁢           (   21   )                 
such that
 
                       T   chn     (     q   sym     )     =           (     -   1     )         (     K   -   1     )     ⁢     q   sym         [         w   chn   *     (         M   ADC     ⁢     q   sym       +     m   ADC       )     ⁢     e       -   j     ⁢   2   ⁢   π   ⁢     k   chn     ⁢       m   ADC     /     M   ADC             ]     .             Eqn   ⁢           (   22   )                 
This channelizer can be implemented using methods well-known to those of ordinary skill in the art.
 
     If K chn =K sub M DoF , then x chn (n sym ) can be directly transformed to K sub  contiguous M DoF ×1 vector subband signals 
               {       x   DoF     (       n   sym     ;     k   sub       )     }         k   sub     =   0         K   sub     -   1           
with contiguous in-subband channels using a M DoF ×K sub  reshaping operation ( 161 ), such that
 
                 x   DoF     (       n   sym     ;     k   sub       )     =         [       x   chn     (           M   DoF     ⁢     k   sub       +     m   DoF       ,     n   sym       )     ]         m   DoF     =   0         M   DoF     -   1       .           
The subband center frequencies
 
               {       f   sub     (     k   sub     )     }         k   sub     =   0         K   sub     -   1           
are then given by
 
                         f   sub     (     k   sub     )     =       (       k   sub     -         K   sub     -   1     2       )     ⁢     M   DoF     ⁢     f   sym         ,   
       k   sub     =   0     ,   …         ,       K   sub     -   1             Eqn   ⁢           (   23   )                 
while the frequency channels within each subband are locally given by
 
                         f   chn     (     m   DoF     )     =       (       m   DoF     -         M   DoF     -   1     2       )     ⁢     f   sym         ,   
       m   DoF     =   0     ,   …         ,       M   DoF     -   1             Eqn   ⁢           (   24   )                 
such that f chn (M DoF k sub +m DoF )=f sub (k sub )+f chn (m DoF ).
 
       FIG.  11    is a block diagram of a single-feed symbol-rate synchronous multi-subband channelizer according to another aspect of the disclosure. In this aspect, the ADC output signal is passed through a set of K sub  frequency converters ( 162   a )-( 162   b ) and decimators ( 163   a )-( 163   b ), each of which frequency downconverts the ADC output signal by frequency f sub (k sub ), and decimates that signal to intermediate rate f dec =M dec f sym , where M dec  is a positive integer. Each decimated signal x dec (n dec ;k sub ) is then passed through a symbol-rate synchronous channelization operation ( 164   a )-( 164   b ) comprising a 1:M dec  serial-to-parallel converter and a M DoF ×M dec  channelization operator T DoF  ( 165 ), to form channelized subband signal x DoF (n sym ;k sub )=T DoF ∘x sym (n sym ;k sub ) where 
                 x   sym     (       n   sym     ;     k   sub       )     =         [       x   dec     (           M   dec     ⁢     n   sym       +     m   dec       ;     k   sub       )     ]         m   dec     =   0         M   dec     -   1       .           
The subband frequencies
 
               {       f   sub     (     k   sub     )     }         k   sub     =   0         K   sub     -   1           
can be organized in accordance with Eqn (23); set to a preselected set of arbitrary frequencies; or adaptively determined based on channel dynamics or co-channel interference considerations.
 
     Given the reception scenario described in  FIG.  1    and  FIG.  2   , and the received signal model given in Eqn (4)-Eqn (14), then the single-sample channelizer output signal shown in  FIG.  7    can be modeled by 
                         x   chn     (     n   sym     )     =         i   chn     (     n   sym     )     +       ∑     ℓ   =   1       L   T           s   chn     (       n   sym     ;   ℓ     )           ,           Eqn   ⁢           (   25   )                 
where i chn (n sym )=T chn ∘i sym (n sym ) and s chn (n sym ;l)=T chn ∘s sym (n sym ;l) are the K chn ×1 interference and navigation signal l channelizer output signals, respectively, and where
 
                   i   sym     (     n   sym     )     =       [       i   R     (       T   sym     (       n   sym     +       m   ADC       M   ADC         )     )     ]     ⁢         and       ⁢   
         s   sym     (       n   sym     ;   ℓ     )     =     [       s   R     (         T   sym     (       n   sym     +       m   ADC       M   ADC         )     ;   ℓ     )     ]             
are the M ADC ×1 interference and GNSS signal l S/P output vectors, respectively.
 
     The &#39;417 NPA, incorporated herein by reference, provides a detailed description of the signals obtaining at the output of a symbol-rate-synchronous channelizer with K chn =M DoF . A simpler channel model can be obtained using the polyphase filter aspect shown in  FIG.  10   . In this case, under the simplifying assumption 
                         s   R     (     t   ;   ℓ     )     ≈         g   TR     (   ℓ   )     ⁢     e     j   ⁢   2   ⁢   π   ⁢       α   TR     (   ℓ   )     ⁢   t       ⁢       s   T     (       t   -       τ   TR     (   ℓ   )       ;   ℓ     )         ,           Eqn   ⁢           (   26   )                 
i.e., ignoring affects of platform dynamics other than Doppler shift caused by the transmitter velocity observed at the receiver, then the channelizer output signal s chn (k chn ,n sym ;l) is given by
 
                         s   chn     (       k   chn     ,       n   sym     ;   ℓ       )     =       e     j   ⁢   2   ⁢   π   ⁢       α   TR     (   ℓ   )     ⁢     T   sym     ⁢     n   sym         ⁢       ∑     q   sym             a   chn     (       k   chn     ,       q   sym     ;   ℓ       )     ⁢       d   T     (         n   sym     -     q   sym       ;   ℓ     )             ,           Eqn   ⁢           (   27   )                 
where a chn (k chn ,q sym ;l) is given by
 
                         a   chn     (       k   chn     ,       q   sym     ;   ℓ       )     =       e     j   ⁢   2   ⁢     π   ⁡   (         q   sym     ⁢     T   sym       -       τ   TR     (   ℓ   )       )     ⁢       f   chn     (     k   chn     )         ⁢     e       -   j     ⁢   2   ⁢   π   ⁢     q   sym     ⁢       α   TR     (   ℓ   )         ×       ∫     -   ∞     ∞           W   chn   *     (     e     j   ⁢   2   ⁢   π   ⁢     fT   ADC         )     ⁢       H   TR     (     f   +       f   chn     (     k   chn     )       )     ⁢     e     j   ⁢   2   ⁢     π   ⁡   (         q   sym     ⁢     T   sym       -       τ   TR     (   ℓ   )       )     ⁢   f       ⁢   df           ,           Eqn   ⁢           (   28   )                 
in which W(e j2πf ) is the discrete Fourier transform of w chn (m ADC ) and H TR (f;l) is the frequency response of the end-to-end symbol shaping, given by
 
 H   TR ( f;l )= g   TR ( l ) j2πα     TR     (l)   H   R ( f ) H   T ( f−α   TR ( l ); l )√{square root over ( P   T ( l ))},  Eqn (29)
 
where H T (f;l) is the Fourier transform of h T (t;l) given in Eqn (2) and H R (f) is the frequency response of the receiver filtering operations.
 
     If H TR (f;l) is bandlimited below f ADC /2 e.g., using typical pre-ADC antialiasing filters, the bandwidth of W(e j2πf ) is much less than the frequency variation in H TR (f;l) and f chn (k chn ) is given by Eqn (21), then Eqn (28) can be approximated by
 
 a   chn ( k   chn   ,q   sym   ;l )≈ a   chn ( k   chn   ;l ) a   sym ( q   sym   ;l )  Eqn (30)
 
where a chn (k chn ;l) and a sym (q sym ;l) are the frequency-varying and time-varying components of the channel signature,
 
                                              a   chn     (       k   chn     ;   ℓ     )     =       f   ADC     ⁢       H   TR     (         f   chn     (     k   chn     )     ;   ℓ     )     ⁢     e       -   j     ⁢   2   ⁢   π   ⁢       f   chn     (     k   chn     )     ⁢       τ   TR     (   ℓ   )             ,             Eqn   ⁢           (   31   )                                   a   sym     (       q   sym     ;   ℓ     )     =         (     -   1     )         (       K   chn     -   1     )     ⁢     q   sym         ⁢       w   sym   *     (       q   sym     -         τ   TR     (   ℓ   )     ⁢     f   sym         )     ⁢     e       -   j     ⁢   2   ⁢   π   ⁢     q   sym     ⁢       α   TR     (   ℓ   )             ,           Eqn   ⁢           (   32   )                 
respectively, referred to herein as the channel frequency signature and channel time signature, and where w sym (τ) is the interpolated channelizer window, given by
 
                       w   sym     (   τ   )     =       ∑     m   ADC             w   chn     (     m   ADC     )     ⁢       Sa   ⁡   (     π   ⁡   (       τ   ⁢     M   ADC       -     m   ADC       )     )     .                 Eqn   ⁢           (   33   )                 
The channelized received signal s chn (k chn ,n sym ;l) is then approximated by
 
 s   chn ( k   chn   ,n   sym   ;l )≈ a   chn ( k   chn   ;l ) d   R ( n   sym   ;l ),  Eqn (34)
 
where d R (n sym ;l) is the signal l symbol sequence observed at the receiver, given by
 
 d   R ( n   sym   ;l )=( a   sym ( n   sym   ;l )∘ d   T ( n   sym   ;l )) e   j2πα     TR     (l)T     sym     n     sym   ,  Eqn (35)
 
and where “∘” here denotes the convolution operation. Further decomposing τ TR (l) and α TR (l) into symbol-normalized TOA and FOA components
 
                       τ   TR     (   ℓ   )     =       (           τ   ~     TR     (   ℓ   )     +       n   TR     (   ℓ   )       )     ⁢     T   sym     ⁢     {                       τ   ~     TR     (   ℓ   )     ∈     [         0       1             )     ,           fractional   ⁢         fine   ⁢         TOA                     n   TR     (   ℓ   )     ∈   ℤ     ,           TOA   ⁢         symbol   ⁢         offset           ,                 Eqn   ⁢           (   36   )                                 α   TR     (   ℓ   )     =       (           α   ~     TR     (   ℓ   )     +       k   TR     (   ℓ   )       )     ⁢     f   sym     ⁢     {                     α   TR     (   ℓ   )     ∈     [           -     1   2             1   2               )     ,           fractional   ⁢         fine   ⁢         FOA                     k   TR     (   ℓ   )     ∈   ℤ     ,           FOA   ⁢         Nyquist   ⁢         zone           ,                 Eqn   ⁢           (   37   )                 
respectively, then Eqn (34) can be replaced by
 
                         d   R     (       n   sym     ;   ℓ     )     =       (         a   sym     (       n   sym     ;   ℓ     )     ∘       d   T     (         n   sym     -       n   TR     (   ℓ   )       ;   ℓ     )       )     ⁢     e     j   ⁢   2   ⁢   π   ⁢         α   ~     TR     (   ℓ   )     ⁢     n   sym             ,           Eqn   ⁢           (   38   )                                   a   sym     (       q   sym     ;   ℓ     )     =         (     -   1     )         (       K   chn     -   1     )     ⁢     q   sym         ⁢       w   sym   *     (       q   sym     -         τ   ~     TR     (   ℓ   )       )     ⁢     e       -   j     ⁢   2   ⁢   π   ⁢     q   sym     ⁢       α   TR     (   ℓ   )             ,           Eqn   ⁢           (   39   )                 
except for phase shift φ(l)=2π(f chn (0)−{tilde over (α)} TR (l))n TR (l), which is subsumed into end-to-end link gain in g TR (l) in Eqn (5). Then d R (n sym ;l) is completely characterized by the channel time signature defined (over fine TOA ranging between 0 and 1; the received TOA measured to integer number of symbols; and the fractional fine FOA ranging between −½ and ½.
 
     For the GNSS signals listed in  FIG.  4   , and assuming |α TR (l)|&lt;&lt;f chp , then a chn (k chn ;l) can be further approximated by 
                         a   chn     (       k   chn     ;   ℓ     )     ≈         M   ADC         M   chp         ⁢     Sa   ⁡   (     π   ⁢       f   chn     (     k   chn     )     ⁢     T   chp       )     ⁢       H   R     (       f   chn     (     k   chn     )     )     ⁢       g   TR     (   ℓ   )     ⁢         P   T     (   ℓ   )       ⁢       ε   R     (       k   chn     ;   ℓ     )         ,           Eqn   ⁢           (   40   )                                   ε   R     (       k   chn     ;   ℓ     )     =       1       M   chp         ⁢       C   T     (       e       -   j     ⁢   2   ⁢     π   ⁡   (         f   chn     (     k   chn     )     -       α   TR     (   ℓ   )       )     ⁢     T   chp         ;   ℓ     )     ⁢     e       -   j     ⁢   2   ⁢     π   ⁡   (           f   chn     (     k   chn     )     -       α   TR     (   ℓ   )       ;   ℓ     )     ⁢       τ   TR     (   ℓ   )             ,           Eqn   ⁢           (   41   )                 
where normalized observed ranging code frequency signature ε R (k chn ;l) captures the cross-channel frequency variability of each channel frequency signature due to the separate ranging code used at each transmitter, and the differing TOA on each transmission path. For well-designed pseudo-random code sequences and M chp &gt;&gt;1, ε R (k chn ;l) can be modeled as a zero-mean, unit-variance, independent and identically-distributed (i.i.d.) complex-Gaussian random process. The channelized output signal can then be approximated by
 
                         x   chn     (     n   sym     )     ≈         i   chn     (     n   sym     )     +       ∑     ℓ   =   1       L   T             a   chn     (   ℓ   )     ⁢       d   R     (       n   sym     ;   ℓ     )             ,           Eqn   ⁢           (   42   )                 
where
 
                 a   chn     (   ℓ   )     =       [       a   chn     (       k   chn     ;   ℓ     )     ]         k   chn     =   0         K   chn     -   1             
is the K chn ×1 signal l channel frequency signature and a chn (k chn ;l) is given by Eqn (40)-Eqn (41) and d R (n sym ;l) is given by Eqn (38).
 
     Assuming the background interference i R (t) given in Eqn (4) is complex-Gaussian and stationary with power spectral density (PSD) S i     R     i     R    (f), then the interference power level in each channel is given by 
                       〈         ❘   &#34;\[LeftBracketingBar]&#34;         i   chn     (       k   chn     ,     n   sym       )       ❘   &#34;\[RightBracketingBar]&#34;       2     〉     =       ∫       -   1     /   2       1   /   2               ❘   &#34;\[LeftBracketingBar]&#34;         W   chn     (     e     j   ⁢   2   ⁢   π   ⁢   f       )       ❘   &#34;\[RightBracketingBar]&#34;       2     ⁢         S   ~         i   R     ⁢     i   R         (     f   +         f   chn     (     k   chn     )     ⁢     T   ADC         )     ⁢   df         ,           Eqn   ⁢           (   43   )                 
where
 
                   S   ~         i   R     ⁢     i   R         (   f   )     =       f   ADC     ⁢       ∑   ℓ             ❘   &#34;\[LeftBracketingBar]&#34;         H   R     (       (     f   +   ℓ     )     ⁢     f   ADC       )       ❘   &#34;\[RightBracketingBar]&#34;       2     ⁢       S       i   R     ⁢     i   R         (       (     f   +   ℓ     )     ⁢     f   ADC       )                 
is the PSD of sampled interference signal i R (T ADC n ADC ). Assuming the ADC input signal is bandlimited to less than f ADC /2 ahead of the sampling operation, and that the PSD frequency variation is low over the passband of the polyphase filter frequency response W chn (e j2πf ), then Eqn (43) can be approximated by
 
                             〈         ❘   &#34;\[LeftBracketingBar]&#34;         i   chn     (       k   chn     ,     n   sym       )       ❘   &#34;\[RightBracketingBar]&#34;       2     〉     ≈           f   ADC     ⁢         ❘   &#34;\[LeftBracketingBar]&#34;         H   R     (       f   chn     (     k   chn     )     )       ❘   &#34;\[RightBracketingBar]&#34;       2     ⁢     S       i   R     ⁢     i   R         ⁢     (       f   chn     ⁢     (     k   chn     )       )     ⁢            ∫         -   1     /   2       1   /   2       ⁢         ❘   &#34;\[LeftBracketingBar]&#34;         W   chn     (     e     j   ⁢   2   ⁢   π   ⁢   f       )       ❘   &#34;\[RightBracketingBar]&#34;       2     ⁢   df                 =           f   ADC     ⁢         ❘   &#34;\[LeftBracketingBar]&#34;         H   R     (       f   chn     (     k   chn     )     )       ❘   &#34;\[RightBracketingBar]&#34;       2     ⁢       S       i   R     ⁢     i   R         (       f   chn     (     k   chn     )     )     ⁢            w   chn          2   2                   =             M   ADC   2       M   chp       ⁢     f   chp     ⁢         ❘   &#34;\[LeftBracketingBar]&#34;         H   R     (       f   chn     (     k   chn     )     )       ❘   &#34;\[RightBracketingBar]&#34;       2     ⁢       S       i   R     ⁢     i   R         (       f   chn     (     k   chn     )     )               ,           Eqn   ⁢           (   44   )                 
where
 
                      w   chn          2   2     =       ∑     m   ADC             ❘   &#34;\[LeftBracketingBar]&#34;         w   chn     (     m   ADC     )       ❘   &#34;\[RightBracketingBar]&#34;       2             
denotes the squared L2 Euclidean norm of the channelizer weights {w chn (m ADC )} ( 158 ), and further assuming normalization constraint ∥w chn ∥ 2   2 =M ADC . More generally, the interference cross-correlation across symbol lag and frequency channel offset can be approximated by
 
                               R       i   chn     ⁢     i   chn         (       k   chn     ,       ℓ   chn     ;     m   sym         )       =   △           〈         i   chn     (       k   chn     ,     n   sym       )     ⁢       i   chn   *     (       ℓ   chn     ,       n   sym     -     m   sym         )       〉                 ≈             M   ADC   2       M   chp       ⁢     f   chp     ⁢       S       i   R     ⁢     i   R         (           f   chn     (     k   chn     )     +       f   chn     (     ℓ   chn     )       2     )                   ×           (     -   1     )         (       K   chn     -   1     )     ⁢     m   sym         ⁢       ρ   chn     (         M   ADC     ⁢     m   sym       ,         k   chn     -     ℓ   chn         M   ADC         )             ,           Eqn   ⁢           (   45   )                 
where
 
                 ρ   chn     (     m   ,   α     )     =       1     M   ADC       ⁢       ∑   n           w   chn     (   n   )     ⁢       w   chn   *     (     n   -   m     )     ⁢     e       -   j     ⁢   2   ⁢   π   ⁢   α   ⁢   n                   
is the discrete-time ambiguity function for channelizer window {w chn (m ADC )} ( 158 ), normalized M ADC  so that |ρ chn (m,α)|≤1.
 
       FIG.  12    depicts the channel time signature ( FIG.  12 A ) and normalized ambiguity function ( FIG.  12 B ) induced by the polyphase-filter based symbol-rate synchronous channelizer ( 160 ) shown in  FIG.  10   , using single-symbol rectangularly-windowed channelizer weights ( 158 ). This Figure also obtains for the FFT-based symbol-rate synchronous channelized shown in  FIG.  8    and  FIG.  9   . As shown in  FIG.  12 A , the channel time signature is only substantive over a single symbol for fine TOA offsets between 0 and 1, and is nearly equal to unity within this region. The received symbol sequence d R (n sym ;l) given in Eqn (38) then simplifies to 
                         d   R     (       n   sym     ;   ℓ     )     ≈         d   T     (         n   sym     -   1   -       n   TR     (   ℓ   )       ;   ℓ     )     ⁢     e     j   ⁢   2   ⁢   π   ⁢         α   ~     TR     (   ℓ   )     ⁢     n   sym             ,           Eqn   ⁢           (   46   )                 
where the phase term −2π{tilde over (α)} TR (l) is subsumed into g TR (l). Similarly, the normalized ambiguity function shown in  FIG.  12 B  is equal to unity at zero channel and symbol offset. Hence the background interference for this channelizer window can be treated as independent between channels and symbols, yielding K chn ×K chn  interference autocorrelation matrix (ACM)
 
                     R       i   chn     ⁢     i   chn         =         M   ADC       M   chp       ⁢     f   chp     ⁢        diag   ⁢              {           ❘   &#34;\[LeftBracketingBar]&#34;         H   R     (       f   chn     (     k   chn     )     )       ❘   &#34;\[RightBracketingBar]&#34;       2     ⁢       S       i   R     ⁢     i   R         (       f   chn     (     k   chn     )     )       }         k   chn     =   0         K   chn     -   1       .               Eqn   ⁢           (   47   )                 
However, the channel model established above can be extended to any channelizer window.
 
     Multiplying x chn (n sym ) by the inverse square-root of R i     chn     i     chn    yields normalized signal 
                                 x   ~     chn     (     n   sym     )     ≈           R       i   chn     ⁢     i   chn           -   1     /   2       ⁢       x   chn     (     n   sym     )         ,                 ≈             ε   chn     (     n   sym     )     +       ∑     ℓ   =   1       L   T             u   chn     (   ℓ   )     ⁢       d   R     (       n   sym     ;   ℓ     )             ,                 Eqn   ⁢           (   48   )                 
where d R (n sym ;l) is given by Eqn (46), ε chn (n sym ) is a K chn ×1 i.i.d. complex-Gaussian random process with zero mean and ACM I K     chn   , and U chn (l) is the K chn ×1 signal l normalized channel frequency signature given by
 
                       u   chn     (   ℓ   )     =       [         Sa   ⁡   (     π   ⁢       f   chn     (     k   chn     )     ⁢     T   chp       )           f   chp     ⁢       S       i   R     ⁢     i   R         (       f   chn     (     k   chn     )     )           ⁢       ε   R     (       k   chn     ;   ℓ     )       ]     ⁢       g   TR     (   ℓ   )     ⁢           P   T     (   ℓ   )       .               Eqn   ⁢           (   49   )                 
In theory, a set of K chn ×1 combining weights {tilde over (w)} max-SINR (l) can then be developed that extracts d R (n sym ;l) from {tilde over (x)} chn (n sym ) with maximum-attainable signal-to-interference-and-noise ratio (max-SINR). For the channel model given in Eqn (48)-Eqn (49), and assuming the navigation symbol sequences are independent for each transmitter, {tilde over (w)} max-SINR (l) is given by
 
                               w     m   ⁢   a   ⁢   x   -   SINR       (   ℓ   )     ∝             (       I     K   chn       +         U   chn     (     ~   ℓ     )     ⁢       U   chn   H     (     ~   ℓ     )         )       -   1       ⁢       u   chn     (   ℓ   )                   =           (       I     K   chn       -         U   chn     (     ~   ℓ     )     ⁢       (       I       L   T     -   1       +       U   chn   H     ⁢     (     ~   ℓ     )     ⁢   
       U   chn     (     ~   ℓ     )         )       -   1       ⁢     U   chn   H         )     ⁢       u   chn     (   ℓ   )               ,           Eqn   ⁢           (   50   )                 
where U chn (˜l)=[u chn (l′)] l′≠l  is the K chn ×(L T −1) matrix of normalized frequency signatures interfering with signal l, and the max-SINR for symbol sequence l is given by
 
                               γ     m   ⁢   a   ⁢   x   -   SINR       (   ℓ   )     ≈                  u   chn     (   ℓ   )          2   2                   -           u   chn   H     (   ℓ   )       ⁢       U   chn     (     ~   ℓ     )     ⁢       (       I       L   T     -   1       +       U   chn   H     ⁢     (     ~   ℓ     )     ⁢       U   chn     
   (     ~   ℓ     )         )       -   1       ⁢       U   chn   H     (     ~   ℓ     )     ⁢       u   chn     (   ℓ   )             .           Eqn   ⁢           (   51   )                 
Modeling ε R (k sub ;l) as i.i.d. complex Gaussian with zero mean and unity variance, appropriate for unknown, well-modeled long ranging codes, then Eqn (51) has mean lower bound
 
                                 γ   _       m   ⁢   a   ⁢   x   -   SINR       (   ℓ   )     =         E   ⁢     {       γ     m   ⁢   a   ⁢   x   -   SINR       (   ℓ   )     }                   ≥         (           γ   _       m   ⁢   a   ⁢   x   -   SNR       (   ℓ   )     -       (       L   T     -   1     )     ⁢     (           γ   _       m   ⁢   a   ⁢   x   -   INR       (   ℓ   )       1   +         γ   _       m   ⁢   a   ⁢   x   -   INR       (   ℓ   )         )     ⁢            max            k   chn       ⁢         γ   _       Rx   -   SNR       (       k   chn     ;   ℓ     )                   ,           Eqn   ⁢           (   52   )                               →   &#34;\[Rule]&#34;       {                   γ   _       m   ⁢   a   ⁢   x   -   SNR       (   ℓ   )     ,                 γ   _       m   ⁢   a   ⁢   x   -   INR       (   ℓ   )     ⁢        1                   (           γ   _       m   ⁢   a   ⁢   x   -   SNR       (   ℓ   )     -       (       L   T     -   1     )            max     k   chn                      γ     _       m   ⁢   a   ⁢   x   -   SNR       ⁢     (       k   chn     ;   ℓ     )         )     ,                   γ   _       m   ⁢   a   ⁢   x   -   INR       (   ℓ   )          ⁢   1           ,               Eqn   ⁢           (   53   )                 
where  γ   Rx-SNR (k chn   ;l) and  γ     max-SNR (l) are the mean receive SNR (Rx-SNR) and maximum-attainable SNR (max-SNR) of signal l in the absence of interference navigation signals, respectively,
 
                                 γ   _       Rx   -   SNR       (       k   chn     ;   ℓ     )     =         E   ⁢     {       ❘   &#34;\[LeftBracketingBar]&#34;         u   chn     (       k   chn     ;   ℓ     )       ❘   &#34;\[RightBracketingBar]&#34;       }                   =               Sa   2     (     π   ⁢       f   chn     (     k   chn     )     ⁢     T   chp       )         f   chp     ⁢       S       i   R     ⁢     i   R         (       f   chn     (     k   chn     )     )         ⁢         ❘   &#34;\[LeftBracketingBar]&#34;         g   TR     (   ℓ   )       ❘   &#34;\[RightBracketingBar]&#34;       2     ⁢       P   T     (   ℓ   )               ,           Eqn   ⁢           (   54   )                                           γ   _       m   ⁢   a   ⁢   x   -   SNR       (   ℓ   )     =         E   ⁢     {              u   chn     (   ℓ   )          2   2     }                   =           ∑       k   chn     =   0         K   chn     -   1             γ   _       Rx   -   SNR       (       k   chn     ;   ℓ     )               ,           Eqn   ⁢           (   55   )                 
and where  γ   max-INR (l) is the mean max-attainable SNR of navigation signals interfering with signal l,
 
                         γ   _       m   ⁢   a   ⁢   x   -   INR       (   ℓ   )     =       1       L   T     -   1       ⁢       ∑       ℓ   ′     ≠   ℓ               γ   _       m   ⁢   a   ⁢   x   -   SNR       (     ℓ   ′     )     .                 Eqn   ⁢           (   56   )                 
As Eqn (53) shows, at high receive SNR, the max-SINR combiner uses L T −1 combiner degrees of freedom (DoF&#39;s) to excise the navigation signals impinging on the receiver, and employs the remaining K chn −L T +1 combiner DoF&#39;s to suppress the background interference i chn (n sym ).
 
     As disclosed in the &#39;417 NPA, incorporated herein by reference, linear algebraic methods to determine max-SINR weights can be devised given either knowledge of the content of {d T (n sym ;l)} l=1   L     T   , e.g., provided by third-party sources; knowledge of a component of {d T (n sym ;l)} l=1   L     T   , e.g., a known/periodic pilot signal; or knowledge of some other exploitable structure of {d T (n sym ;l)} l=1   L     T   , as listed in  FIG.  4   . These methods can be used to detect the navigation signals in the environment, estimate the max-SINR of the combiner output symbol sequences for those signals, and determine at least the geo-observables {n TR (l),{tilde over (α)} TR (l)} l=1   L     T    for those signals, and if needed estimate information-bearing components of the transmitted symbol streams {d T (n sym ;l)} l=1   L     T   . Moreover, these techniques can perform these operations without prior knowledge of the signal ranging code or receive frequency signature for the signals, or a search over TOA and FOA is the ranging code is known. 
     However, as also disclosed in the &#39;417 NPA, linear-algebraic methods that can develop combiner weights covering the entire navigation signal passband, e.g., K chn  ˜1,000 for the GPS L1 legacy signal, and K chn  ˜10,000 for the GPS L5 civil signal, require O(K chn   2 ) operations/symbol to be implemented. Moreover, they require O(K chn ) symbols to converge to a useful solution, e.g., 2-to-4 seconds for the GPS L1 legacy signal, and 20-40 seconds for the GPS L5 civil signal. This convergence time introduces additional complexity in GNSS reception scenarios, due to channel dynamics caused by movement of the GNSS satellite vehicles. 
     In the &#39;417 NPA, partial subband methods that reduce K chn  to a manageable value are disclosed to overcome these issues. These methods can provide strong advantage in the presence of strong MOS-DSSS spoofing and jamming signals. However, they sacrifice substantive processing gain to achieve this capability at reasonable complexity and over short reception intervals. This issue can be especially significant for modern wideband navigation signals such as the GPS L5 civil signal. The multi-subband approach provides a path to overcome this limitation. 
     Returning to  FIG.  10   , reshape operation ( 161 ) shown in  FIG.  10    separates x chn (n sym ) into K sub  subbands, each comprising M DoF  frequency channels such that K chn =K sub M DoF , and the signal in each subband can be represented by 
                               x   DoF     (       n   sym     ;     k   sub       )     =         (             x   chn     (         M   DoF     ⁢     k   sub       ,     n   sym       )             ⋮               x   chn     (           M   DoF     ⁢     k   sub       +     M   DoF     -   1     ,     n   sym       )           )                 ≈             i   DoF     (       n   sym     ;     k   sub       )     +       ∑     ℓ   =   1       L   T             a   DoF     (       k   sub     ;   ℓ     )     ⁢       d   R     (       n   sym     ;   ℓ     )                   ,           Eqn   ⁢           (   57   )                 
where d R (n sym ;l), given by Eqn (46) for the rectangularly-windowed channelizer, is shared for all subbands, and a DoF (k sub ;l) is the M DoF ×1 signal l frequency signature for subband k sub ,
 
                         a   DoF     (       k   sub     ;   ℓ     )     =     (             a   chn     (         M   DoF     ⁢     k   sub       ,   ℓ     )             ⋮               a   chn     (           M   DoF     ⁢     k   sub       +     M   DoF     -   1     ,   ℓ     )           )       ,           Eqn   ⁢           (   58   )                 
and where i DoF (n sym ;k sub ) is the background noise and interference received in subband k sub . If i R (t) is stationary with PSD S i     R     i     R   (f) at the LPF output, then the zero-lag ACM of i DoF (n sym ;k sub ) is further approximated by
 
                         R       i   DoF     ⁢     i   DoF         (     k   sub     )     ≈         M   ADC   2       M   chp       ⁢     f   chp     ⁢         ❘   &#34;\[LeftBracketingBar]&#34;         H   R     (       f   sub     (     k   sub     )     )       ❘   &#34;\[RightBracketingBar]&#34;       2     ⁢       S       i   R     ⁢     i   R         (       f   sub     (     k   sub     )     )     ⁢     P   sub         ,           Eqn   ⁢           (   59   )                 
where
 
                 f   sub     (     k   sub     )     =       (       k   sub     -         K   sub     -   1     2       )     ⁢     M   DoF     ⁢     f   sym     ⁢         i           
is the center frequency of subband k sub , and where
 
               P   sub     =       [       1     M   ADC       ⁢       ∑   m             ❘   &#34;\[LeftBracketingBar]&#34;         w   chn     (   m   )       ❘   &#34;\[RightBracketingBar]&#34;       2     ⁢     e       -   j     ⁢   2   ⁢     π   ⁡   (     k   -   ℓ     )     ⁢   m   /     M   ADC               ]       k   ,     ℓ   =   0           M   DoF     -   1             
is the normalized subband interference ACM, which is equal to I M     DoF    for the single-symbol rectangular channelizer window.
 
     If M DoF ≥L T , then navigation symbol sequences l can be received in each subband using an appropriately designed set of max-SINR linear in-subband linear combining weights w max-SINR (k sub ;l), such that combiner output symbol sequence {circumflex over (d)} R (n sym ;k sub ;l)=w max-SINR   H (k sub ;l)x sub (n sym ;k sub ) excises the interfering navigation signals {d R (n sym ;l′)} l′≠l  also received in the subband using L T −1 degrees of freedom, and suppresses the background interference i sub (n sym ;k chn ) in the subband using the remaining M DoF −L T +1 degrees of freedom, and where (⋅) H  denotes the Hermitian transpose operation. Assuming a nearly-flat signature and interference response over each subband, and rectangular channelizer windows, then the SINR achievable using this linear combiner satisfies mean lower bound 
                           γ   _       m   ⁢   a   ⁢   x   -   SINR       (       k   sub     ;   ℓ     )     ≥       (       M   DoF     -       (       L   T     -   1     )     ⁢     (           γ   _       m   ⁢   a   ⁢   x   -   INR       (       k   sub     ;   ℓ     )       1   +         γ   _       m   ⁢   a   ⁢   x   -   INR       (       k   sub     ;   ℓ     )         )         )     ⁢         γ   _       Rx   -   SNR       (       k   sub     ;   ℓ     )         ,           Eqn   ⁢           (   60   )                                                    →   &#34;\[Rule]&#34;       {               M   DoF     ⁢         γ   _       Rx   -   SNR       (       k   sub     ;   ℓ     )       ,                 γ   _       m   ⁢   a   ⁢   x   -   INR       (       k   sub     ;   ℓ     )     ⁢        1                     (       M   DoF     -     L   T     +   1     )     ⁢         γ   _       Rx   -   SNR       (       k   sub     ;   ℓ     )       ,                   γ   _       m   ⁢   a   ⁢   x   -   INR       (       k   sub     ;   ℓ     )          ⁢   1                       Eqn   ⁢           (   61   )                 
in each subband, where  γ   Rx-SNR (k sub ;l) and  γ   max-INR (k sub ;l) are the mean receive signal l SNR and the max-attainable navigation signal interference received along with that signal in subband k sub , respectively
 
                           γ   _       Rx   -   SNR       (       k   sub     ;   ℓ     )     ≈           Sa   2     (     π   ⁢       f   sub     (     k   sub     )     ⁢     T   chp       )         f   chp     ⁢       S       i   R     ⁢     i   R         (       f   sub     (     k   sub     )     )         ⁢         ❘   &#34;\[LeftBracketingBar]&#34;         g   TR     (   ℓ   )       ❘   &#34;\[RightBracketingBar]&#34;       2     ⁢       P   T     (   ℓ   )         ,           Eqn   ⁢           (   62   )                                   γ   _       m   ⁢   a   ⁢   x   -   INR       (       k   sub     ;   ℓ     )     ≈         M   DoF         L   T     -   1       ⁢       ∑       ℓ   ′     ≠   ℓ               γ   _       Rx   -   SNR       (       k   sub     ;     ℓ   ′       )     .                 Eqn   ⁢           (   63   )                 
As Eqn (61) shows, the max-SINR weights can excise up to L T −1 strong MOS-DSSS signals impinging on the receiver in each subband, and can provide an additional factor of M DoF −L T +1 SNR gain to suppress the remaining background noise and interference in each subband. In low-SNR environments where all of the navigation signals are received well below the noise floor ( γ   RX-SNR (k sub ;l)&lt;&lt;1/M DoF ), the max-SINR weights can provide the full factor-of-M DoF  SNR gain available to the receiver.
 
     In aspects of the invention, the symbol sequence estimates from each subband are then further combined using cross-band linear combining weights 
                 {       g   sub     (       k   sub     ;   ℓ     )     }         k   sub     =   0         K   sub     -   1       ,         
yielding multi-subband symbol sequence estimate
 
                                 d   ^     R     (       n   sym     ;   ℓ     )     ≈           ∑       k   sub     =   0         K   sub     -   1             g   sub   *     (       k   sub     ;   ℓ     )     ⁢         d   ^     R     (       n   sym     ;     k   sub     ;   ℓ     )                     ≈             g   sub   H     (   ℓ   )     ⁢         d   ^     sub     (       n   sym     ;   ℓ     )               .           Eqn   ⁢           (   64   )                 
Assuming each subband has substantively excised the MOS-DSSS signals impinging on the receiver, and has provided a unit-power output signal, then {circumflex over (d)} sub (n sym ;l) can be approximated by
 
                           d   ^     R     (       n   sym     ;   ℓ     )     ≈             ∑     m   ⁢   a   ⁢   x   -   SINR       1   /   2           (   ℓ   )     ⁢     ε   ⁡   (       n   sym     ;   ℓ     )         +         a     m   ⁢   a   ⁢   x   -   SINR       (   ℓ   )     ⁢         d   ^     R     (       n   sym     ;   ℓ     )                         ∑     m   ⁢   a   ⁢   x   -   SINR         (   ℓ   )       =         diag   ⁢     {         γ     m   ⁢   a   ⁢   x   -   SINR       (       k   sub     ;   ℓ     )         (     1   +       γ     m   ⁢   a   ⁢   x   -   SINR       (       k   sub     ;   ℓ     )       )     2       }                       a     m   ⁢   a   ⁢   x   -   SINR       (   ℓ   )     =         [         γ     m   ⁢   a   ⁢   x   -   SINR       (       k   sub     ;   ℓ     )       1   +       γ     m   ⁢   a   ⁢   x   -   SINR       (       k   sub     ;   ℓ     )         ]             .         
The SINR of this signal is maximized by setting g sub (l)=[1+γ Max-SINR (k sub ;l)], yielding output SINR
 
                         γ     m   ⁢   a   ⁢   x   -   SINR       (       n   sym     ;   ℓ     )     =       ∑       k   sub     =   0         K   sub     -   1           γ     m   ⁢   a   ⁢   x   -   SINR       (       k   sub     ;   ℓ     )         ,           Eqn   ⁢           (   65   )                 
which has mean lower bound
 
                           γ   _       m   ⁢   a   ⁢   x   -   SINR       (   ℓ   )     ≥       ∑       k   sub     =   0         K   sub     -   1           (       M   DoF     -       (       L   T     -   1     )     ⁢     (           γ   _       m   ⁢   a   ⁢   x   -   INR       (       k   sub     ;   ℓ     )       1   +         γ   _       m   ⁢   a   ⁢   x   -   INR       (       k   sub     ;   ℓ     )         )         )     ⁢         γ   _       Rx   -   SNR       (       k   sub     ;   ℓ     )           ,           Eqn   ⁢           (   66   )                               →   &#34;\[Rule]&#34;       {                 M   DoF     ⁢       ∑     k   sub             γ   _       Rx   -   SNR       (       k   sub     ;   ℓ     )         ,                 γ   _       m   ⁢   a   ⁢   x   -   INR       (       k   sub     ;   ℓ     )     ⁢        1                     (       M   DoF     -     L   T     +   1     )     ⁢       ∑     k   sub             γ   _       Rx   -   SNR       (       k   sub     ;   ℓ     )         ,                   γ   _       m   ⁢   a   ⁢   x   -   INR       (       k   sub     ;   ℓ     )          ⁢   1           .               Eqn   ⁢           (   67   )                 
Thus, the multi-subband solution can recover all or most of the processing gain of the full-channelizer max-SINR combiner, if M DoF  is greater than the expected number of MOS-DSSS signals impinging on the receiver.
 
     In-subband linear combining weights approaching the max-SINR solution can be computed using linear-algebraic methods disclosed in the &#39;417 NPA, given either knowledge of the content of {d T (n sym ;l)} l=1   L     T   , e.g., provided by third-party sources; knowledge of a component of {d T (n sym ;l)} l=1   L     T   , e.g., a known/periodic pilot signal; or other exploitable structure of {d T (n sym ;l)} l=1   L     T   , as listed in  FIG.  4   . These methods also can be used to estimate the max-SINR obtaining in each subband, thereby providing cross-subband weights for the second combining stage shown in Eqn (64). The multi-subband solution can further be used to estimate the full max-SINR of the received symbol sequences, detect the navigation signals in the environment, determine at least the geo-observables {n TR (l),{tilde over (α)} TR l)} l=1   L     T    for those signals, and if needed estimate information-bearing components of the transmitted symbol streams {d T (n sym ;l)} l=1   L     T   . Moreover, these techniques can perform these operations without prior knowledge of the signal ranging code or receive frequency signature for the signals, or a search over TOA and FOA is the ranging code is known. 
     In contrast to the full-band solution, however, the complexity of the multi-subband method scales quadratically with M DoF  and linearly with K sub , such that the full set of K chn  channels provided by the symbol-rate synchronous channelizer can be processing at a complexity of O(K sub M DoF   2 )=O(K chn M DoF ) operations/symbol—a factor of K sub  saving in operations relative to the full-band system. Moreover, the number of symbols needed to provide a stable solution using the methods disclosed in the &#39;417 NPA are O(M DoF ), also a factor of K sub  reduction in convergence time relative to the full-band methods. In the invention, the value of M DoF  becomes a design parameter dictated by the number of MOS-DSSS signals expected to be encountered in the receive environment, and the desired complexity and convergence time of the signal reception algorithms, independent from the bandwidth of the receiver. 
     As an specific example, setting K chn =960 for the GPS L1 legacy signal and K chn =9,600 for the GPS L5 civil signal, i.e., and setting M DoF =60 in both cases to allow reception of all of the GPS signals in the field of view of the receiver as well as 30-45 beacons or spoofers, then the multi-subband processor has a complexity and convergence time reduction of 16 for the GPS L1 legacy signal, and 160 for the GPS L5 civil signal, allowing convergence in as little as 120 ms for both methods and complexity of O(57.6) Mops/s for the GPS L1 legacy signal and O(576) Mops/s for the GPS L5 civil signal—well within the capabilities of modern DSP equipment, and low-cost DSP for the GPS L1 legacy signal. Moreover, the multi-subband approach is highly amenable to parallel processing methods, facilitating it&#39;s implementation in FPGA&#39;s and general-purpose GPU (GPGPU) architectures. 
     Taking channel time-variation given in Eqn (5)-Eqn (14) due to moving transmitters and/or receivers more precisely into account, then over reception intervals on the order of 0-to-2 seconds, Eqn (42) generalizes to 
                         x   chn     (     n   sym     )     ≈         i   chn     (     n   sym     )     +       ∑     ℓ   =   1       L   T             a   chn     (       n   sym     ;   ℓ     )     ⁢       d   R     (       n   sym     ;   ℓ     )             ,           Eqn   ⁢           (   68   )                 
where time-varying frequency-signature a chn (n sym ;l) is approximated by
 
                         a   chn     (       n   sym     ;   ℓ     )     ≈         a   chn     (   ℓ   )     ⊙     [     e       -   j     ⁢   2   ⁢     π   ⁡   (         f   chn     (     k   chn     )     -       α   TR     (   ℓ   )       )     ⁢       T   sym     (           τ   ~     TR     (   ℓ   )     +         τ   TR     (   1   )       (   ℓ   )     ⁢     n   sym         )         ]         ,           Eqn   ⁢           (   69   )                 
and the observed symbol sequence is given by
 
                               d   R     (       n   sym     ;   ℓ     )     ≈               (       a   sym     ⁢       (       n   sym     ;   ℓ     )     ⊗                         d   T     ⁢     (         n   sym     -       n   TR     ⁢     (   ℓ   )         ;   ℓ     )       )     ⁢     e     j   ⁢   2   ⁢     π   ⁡   (             α   ~     TR     (   ℓ   )     ⁢     n   sym       +       1   2     ⁢         α   ~     TR     (   1   )       (   ℓ   )     ⁢     n   sym   2         )                             →             d   T     (         n   sym     -       n   TR     (   ℓ   )     -   1     ;   ℓ     )     ⁢     e     j   ⁢   2   ⁢     π   ⁡   (             α   ~     TR     (   ℓ   )     ⁢     n   sym       +       1   2     ⁢         α   ~     TR     (   1   )       (   ℓ   )     ⁢     n   sym   2         )                   ,           Eqn   ⁢           (   70   )                 
and where “⊙” denotes the element-wise multiplication (Hadamard product) operation, τ TR   (1) (l) is the observed differential TOA (DTOA) given in Eqn (10), and {acute over (α)} TR   (1) (l)=α TR   (1) (l)T sym   2  is the symbol-normalized observed differential FOA (DFOA). Hence the DTOA principally affects the channel frequency signature, by inducing a linear frequency shift across the symbol-rate synchronous frequency channels, while the DFOA principally affects the channel time signature, by inducing a quadratic frequency shift over the time symbols. However, it should be noted that the FOA is linearly related to the DTOA, as shown in Eqn (9), hence it affects both signature components.
 
     For typical GNSS reception scenarios, τ TR   (1) (l) is typically less than 4 μs/s in magnitude, resulting in a link FOA of ±6.3 kHz at the L1 GPS center frequency, and ±4.7 kHz at the GPS L5 center frequency. From Eqn (69), This DTOA value can also cause a differential frequency shift of as much as 8 Hz over the ±1.023 MHz L1 GPS legacy signal band, or 80 Hz over the ±10.23 MHz L5 civil GPS band; hence it must be taken into account in any geo-observable estimation algorithm implemented over &gt;10 ms reception intervals. 
     In this regard, the multi-subband approach shown in  FIG.  10    provides an additional benefit, by separating the frequency channels into narrow subbands during the reshape operation ( 161 ). In this case, the subband signal can be modeled by 
                                              x   DoF     (     n   sym     )     ≈         i   DoF     (     n   sym     )     +       ∑     ℓ   =   1       L   T             a   DoF     (       n   sym     ;   ℓ     )     ⁢       d   R     (       n   sym     ;     k   sub     ;   ℓ     )             ,             Eqn   ⁢           (   71   )                                 d   R     (       n   sym     ;     k   sub     ;   ℓ     )     =           Eqn   ⁢           (   72   )                             d   T     (         n   sym     -       n   TR     (   ℓ   )     -   1     ;   ℓ     )     ⁢     e       -   j     ⁢   2   ⁢     π   ⁡   (       f   T     +       f   sub     (     k   sub     )       )     ⁢       T   sym     (           τ   TR     (   1   )       (   ℓ   )     ⁢     n   sym       +       1   2     ⁢         τ   ~     TR     (   2   )       (   ℓ   )     ⁢     n   sym   2         )           ,                         a   DoF     (       n   sym     ;   ℓ     )     =           a   DoF     (   ℓ   )     ⊙     [     e       -   j     ⁢   2   ⁢     π   ⁡   (       m   DoF     -     (         M   DoF     -   1     2     )     -         α   TR     (   ℓ   )     ⁢     T   sym         )     ⁢     (           τ   ~     TR     (   ℓ   )     +         τ   TR     (   1   )       (   ℓ   )     ⁢     n   sym         )         ]       .             Eqn   ⁢           (   73   )                 
In aspects of the invention described herein, the substantive cross-subband frequency shift modeled in Eqn (72) is estimated as part of the geo-observable estimation procedure. The observed time-variability in the subband frequency signature a DoF (n sym ;l) modeled in Eqn (73) induces a much lower in-subband frequency shift, due to the much lower bandwidth of the subband. For the example cited above, where M DoF =60, a DTOA of 4 μs/s induces a frequency shift of 0.24 Hz across the subband, or 0.12 cycles over 500 ms. As disclosed in the &#39;417 NPA, this frequency shift creates low-level dispersive components that can be excised as part of the adaptation process for sufficiently large values of M DoF ; the value considered here is likely to be more than sufficient for typical GNSS reception scenarios, even in the presence of strong ground beacons, and in the presence of strong spoofers attempting to emulate GNSS transmitter dynamics.
 
     Similarly, for typical GNSS reception scenarios, α TR   (1) (l) is typically less than 3 Hz/s in magnitude. As a consequence, its affect is minimal for reception intervals of 500 ms or less, and manageable for reception intervals on the order of 2 seconds or less, using aspects of the invention disclosed herein. 
     Accounting for worst-case channel dispersion, if M DoF ≥2L T  the entire network symbol stream can be extracted from the channel, i.e., the interfering MOS-DSSS symbol sequences can be excised from each GNSS symbol, using purely linear combining operations, and even if the navigation signals are received with high-SNR inducing significant cross-interference. 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&#39;s, and in the absence of pseudolites or spoofers operating in the same band as those SV&#39;s, no more than 12 SV&#39;s are likely to be within the field of view of a GNSS receiver at any one time (L T ≤12). Over short observation intervals, maximum substantive rank change induced by the MOS-DSSS signals is therefore likely to be 24. This number can grow to 48 in the presence of 12 spoofers “assigned” to each legitimate navigation signal. In both cases, this is much less than the number of channels available for any GNSS signals listed in  FIG.  4   . 
     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.  4   , 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 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). 
     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&#39;s and Kryder&#39;s Laws, the ability to fully exploit this channel response will become increasingly easier over time. 
     The previous description is provided to enable any person skilled in the art to practice the various aspects described herein. Various modifications to these aspects will be readily apparent to those skilled in the art, and the generic principles defined herein may be applied to other aspects. Thus, the claims are not intended to be limited to the aspects shown herein, but is to be accorded the full scope consistent with the language claims, wherein reference to an element in the singular is not intended to mean “one and only one” unless specifically so stated, but rather “one or more.” The word “exemplary” is used herein to mean “serving as an example, instance, or illustration.” Any aspect described herein as “exemplary” is not necessarily to be construed as preferred or advantageous over other aspects. Unless specifically stated otherwise, the term “some” refers to one or more. The phrase “A or B” may correspond to A only, B only, or A and B. Combinations such as “at least one of A, B, or C,” “one or more of A, B, or C,” “at least one of A, B, and C,” “one or more of A, B, and C,” and “A, B, C, or any combination thereof” include any combination of A, B, and/or C, and may include multiples of A, multiples of B, or multiples of C. Specifically, combinations such as “at least one of A, B, or C,” “one or more of A, B, or C,” “at least one of A, B, and C,” “one or more of A, B, and C,” and “A, B, C, or any combination thereof” may be A only, B only, C only, A and B, A and C, B and C, or A and B and C, where any such combinations may contain one or more member or members of A, B, or C. The words “module,” “block,” “element,” “device,” and the like may not be a substitute for the word “means.” As such, no claim element is to be construed as a means plus function unless the element is expressly recited using the phrase “means for.” 
     Extension to Multifeed Reception 
       FIG.  13    depicts aspects of the invention, in which the MOS-DSSS navigation signals are received over M feed  feeds using a set of spatial and/or polarization diverse antennas ( 90   a - 90   c ). In some aspects, that antennas are additionally coupled from a larger set of antennas to M feed  feeds using a beamforming network (BFN) ( 94 ), e.g., a Butler mode former, to provide spatial beamsteering operations, exclude known interference sources, or extract stronger and/or spatially whitened modes from the antenna array. The M feed  feeds are then coherently converted to complex-baseband using a direct-conversion receiver bank ( 200 ), to form M feed ×1 signal vector 
                 x   R     (   t   )     =         [       x   R     (     t   ;     m   feed       )     ]         m   feed     =   1       M   feed       .           
The feed vector x R (t) is then sampled by a bank of dual-ADC&#39;s ( 243 ), at sampling rate M ADC f sym  determined by a common clock ( 117 ) where f sym  is the baseband symbol rate of the MOS-DSSS navigation signals transmitted to the receiver, and M ADC  is a positive integer.
 
     Assuming coherent downconversion of the received signals, the complex M feed ×1 signal generated output from the receiver bank ( 200 ) is modeled as 
                   x   R     (   t   )     =         i   R     (   t   )     +       ∑     ℓ   =   1     L             s   R     (     t   ;   ℓ     )           ,         
where i R (t) comprises the M feed ×1 vector of background noise and co-channel interference (CCI) present in the receiver passband, and s R (t;l) is M feed ×1 vector of MOS-DSSS navigation signals received from transmitted l. In the absence of nonlocal multipath, S R (t;l) is modeled by
 
                         s   R     (         t   0     +   t     ;   ℓ     )     ≈       e     j   ⁢   2   ⁢     π   ⁡   (           α   TR     (   ℓ   )     ⁢   t     +       1   2     ⁢       α   TR     (   1   )       (   ℓ   )     ⁢     t   2         )         ⁢       g   TR     (   ℓ   )     ⁢     (         a   TR     (     t   ;   ℓ     )     ∘       s   T     (           (     1   -       τ   TR     (   1   )       (   ℓ   )       )     ⁢   t     -       τ   TR     (   ℓ   )       ;   ℓ     )       )         ,           Eqn   ⁢           (   74   )                 
where a TR (t;l) is the M feed ×1 time-varying spatial signature operator with frequency response A TR (f;l), which is also assumed to adhere to a first-order spatial signature blur model in presence of channel dynamics, due to both movement of the receiver over the reception interval, and adjustments in the yaw, pitch, and roll orientation of the receiver platform. Assuming the antennas have nonzero gain along right-hand and left-hand circular polarizations, and in absence of local scattering multipath and channel dynamics, a TR (t;l) can be modeled as
 
 A   TR ( f;l )= A   R ( f;Ψ   R   u   TR ( l ))ρ TR ( l ),  Eqn (75)
 
where u TR (l) is the observed line-of-bearing (LOB) direction vector given in Eqn (15) and Ψ R  is a 3×3 rotation matrix that captures the yaw, pitch, and tilt of the receiver platform and converts the LOB vector to a local unit-norm direction-of-arrival (DOA) vector, in some aspects parameterized with respect to azimuth and elevation angles, and where
 
                           {       A   R     (     f   ;   u     )     }              u        2     =   1       =     {     [         A   R     (     f   ;   u     )         ❘   &#34;\[LeftBracketingBar]&#34;     RHCP                 A   R     (     f   ;   u     )               ❘   &#34;\[RightBracketingBar]&#34;       LHCP     ]     }              u        2     =   1           
is the array manifold of M 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 ρ TR (l) is a 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.
 
     As shown in  FIG.  13   , each element of signal vector x R (t) is then channelized into a K chn ×1 signal vector 
                   x   chn     (       n   sym     ;     m   feed       )     =       [       x   chn     (       k   chn     ,       n   sym     ;     m   feed         )     ]         k   chn     =   0         K   chn     -   1         ,         
using a bank of M feed  symbol-rate synchronous 1:K chn  vector channelizers ( 250 ), for example, using M feed  1:M ADC  S/P operations and K chn ×M ADC  channelization matrix operators T chn  ( 253 ), to form channelized signal x chn (n sym ;m feed )=T chn ∘[x R (n sym T sym +m ADC T ADC ;m feed )]. The M feed  channelizer output signals
 
               {       x   chn     (       n   sym     ;     m   feed       )     }         m   feed     =   1       M   feed           
are then combined to form a single M chn ×1 complex vector x chn (n sym ) ( 259 ), where M chn =K chn M feed . In some aspects of the invention, this is performed by “stacking” the signals first by receiver feed, and then by channel, such that
 
                       x   chn     (     n   sym     )     =       (             [       x   chn     (     0   ,       n   sym     ;     m   feed         )     ]         m   feed     =   1       M   feed               ⋮               [       x   chn     (         K   chn     -   1     ,       n   sym     ;     m   feed         )     ]         m   feed     =   1       M   feed             )     .             Eqn   ⁢           (   76   )                 
In others, combining is performed first by channel, then by receiver feed, such that
 
                 x   chn     (     n   sym     )     =         [       x   chn     (       n   sym     ;     m   feed       )     ]         m   feed     =   1       M   feed       .           
Other, more general combining strategies can also be used; however, the form shown in Eqn (76) is has strong advantages for multi-subband processing.
 
       FIG.  14    depicts one aspect for performing multifeed, multi-subband channelization, using an FFT-based frequency channelization method. The M feed ×1 vector ADC output signal x ADC (n ADC )=x R (T ADC n ADC ) is passed through a bank of M feed  1:M ADC  serial-to-parallel (S/P) convertors ( 251 ) and optionally windowed ( 153 ) symbol-rate synchronous FFT operations ( 252 ), each of which computes a K FFT −bin FFT using M ADC  contiguous samples of ADC output data, where f ADC =M ADC f sym . The FFT operation creates K FFT M feed ×1 data vector 
               [       x   FFT     (       k   bin     ,     n   sym       )     ]         k   bin     =   0         K   FFT     -   1           
with rate f sym  where
 
     
       
         
           
             
               
                 x 
                 FFT 
               
               ( 
               
                 
                   k 
                   bin 
                 
                 , 
                 
                   n 
                   sym 
                 
               
               ) 
             
             = 
             
               
                 
                   [ 
                   
                     
                       x 
                       FFT 
                     
                     ( 
                     
                       
                         k 
                         bin 
                       
                       , 
                       
                         
                           n 
                           sym 
                         
                         ; 
                         
                           m 
                           feed 
                         
                       
                     
                     ) 
                   
                   ] 
                 
                 
                   
                     m 
                     feed 
                   
                   = 
                   0 
                 
                 
                   M 
                   feed 
                 
               
               . 
             
           
         
       
     
     The FFT output vector is then passed through a subband selection operation ( 254 ), which selects K sub  subbands of FFT bins to form subbands 
                 {       x   DoF     (       n   sym     ;     k   sub       )     }         k   sub     =   0         K   sub     -   1       ,         
each with rate f sym . In the aspect shown in  FIG.  14   , each of the subbands comprise M sub  FFT bins, such that x DoF (n sym ;k sub ) is an
 
                   M   DoF     ×   1   ⁢         vector   ⁢         given   ⁢         by   ⁢             x   DoF     (       n   sym     ;     k   sub       )       =       [       x   FFT     (         k   bin     (       m   sub     ;     k   sub       )     ,     n   sym       )     ]         m   sub     =   0         M   sub     -   1         ,         
and where M DoF =M sub M feed . Other aspects may vary the number of FFT bins in each subband, e.g., to account for interference loading or spectral shaping issues in each subband. In the aspect shown in  FIG.  14   , a “dense” subband selection strategy in which k bin (m sub ;k sub )=(k init (k sub )+m sub )mod K FFT  is also shown; however, as in the single-feed receiver, different strategies can be used to form each subband set.
 
       FIG.  15    depicts another aspect of the invention, in which the frequency channelization is performed using a bank of polyphase filter based symbol-rate synchronous channelizers. The M feed ×1 vector ADC output signal x ADC (n ADC )=x R (T ADC n ADC ) is passed through a bank of M feed  1:M ADC  serial-to-parallel (S/P) convertors ( 251 ) and a bank of M feed  polyphase filter based channelizers ( 255 ), each using channelizer weights {w chn (m ADC )} ( 158 ), which creates a K chn M feed ×1 data vector 
               [       x   chn     (       k   chn     ,     n   sym       )     ]         k   chn     =   0         K   chn     -   1           
with rate f sym , where
 
                 x   chn     (       k   chn     ,     n   sym       )     =         [       x   chn     (       k   chn     ,       n   sym     ;     m   feed         )     ]         m   feed     =   0       M   feed       .           
As in the single-feed channelizer shown in  FIG.  10   , the frequency channels are given by Eqn (21).
 
     The K chn M feed ×1 data vector is then put through a M DoF ×K sub  reshape operation ( 161 ), where K chn =K sub M sub , and M DoF =M sub M feed , such that K chn M feed =K sub M DoF . This creates K sub  subband signals 
                 {       x   DoF     (       n   sym     ;     k   sub       )     }         k   sub     =   0         K   sub     -   1       ,         
where M DoF ×1 subband k sub  signal
 
     
       
         
           
             
               
                 x 
                 DoF 
               
               ( 
               
                 
                   n 
                   sym 
                 
                 ; 
                 
                   k 
                   sub 
                 
               
               ) 
             
             = 
             
               
                 
                   [ 
                   
                     
                       x 
                       chn 
                     
                     ( 
                     
                       
                         
                           
                             K 
                             sub 
                           
                           ⁢ 
                           
                             k 
                             sub 
                           
                         
                         + 
                         
                           m 
                           sub 
                         
                       
                       , 
                       
                         n 
                         sym 
                       
                     
                     ) 
                   
                   ] 
                 
                 
                   
                     m 
                     sub 
                   
                   = 
                   0 
                 
                 
                   
                     M 
                     sub 
                   
                   - 
                   1 
                 
               
               . 
             
           
         
       
     
       FIG.  16    and  FIG.  17    show exemplary values of K chn  and K sub , and M sub , designed to yield a constant value of M DoF =60 as the number of receiver feeds is varied from 1 to 6, for the GPS L1 legacy signal ( FIG.  16   ) and the GPS L5 civil signal ( FIG.  17   ). In all cases, the receiver channelizer spans roughly 96% of the signal mainlobe, or just past the half-power signal bandwidth. As the number of receiver feeds grows, the number of frequency channels per subband is dropped, and the number of subbands is increased correspondingly, so that the degrees of freedom for linear-algebraic processing is the same in each subband, and the multi-subband channelizer continues to occupy the same overall bandwidth. 
     The order-of-magnitude complexity in Mop/second and convergence time in milliseconds is also shown in  FIG.  16    and  FIG.  17   . As these Figures show, the receiver complexity grows linearly with the number of receiver feeds and bandwidth, and is well within the capabilities of low-cost DSP for the L1 legacy signal, and moderate cost DSP for the GPS L5 signal. In contrast, the needed convergence time can be a low multiple of 60 ms for both signals and all receiver feeds. Moreover, the value of  FIG.  16    and  FIG.  17    chosen here should be sufficient to detect and demodulate all of the GPS signals within the field of view of any terrestrial receiver, with margin left over for detection and removal of pseudolites and spoofers possibly impinging on the GPS band. Lastly, the multifeed receiver should be able excise as many as M feed −1 wideband non-GPS signals, e.g., jammers, or as many as K sub −1 narrowband jammers. 
     The channel model developed for the single-feed receiver also extends in a straightforward fashion to the multifeed receiver. Ignoring channel dynamics, the polyphase-filter based channelizer yields the same channelized receive signal model given in Eqn (25), and the subband model given in Eqn (58), where 
                   	           a   chn     (   ℓ   )     ≈     [         a   chn     (         f   chn     (     k   chn     )     ;   ℓ     )     ⁢       A   TR     (         f   chn     (     k   chn     )     ;   ℓ     )       ]       ,             Eqn   ⁢           (   77   )                                   a   DoF     (       k   sub     ;   ℓ     )     ≈         [       a   chn     (         f   chn     (       K   sub     +     m   sub       )     ;   ℓ     )     ]         m   sub     =   0         M   sub     -   1       ⊗       A   TR     (         f   sub     (     k   sub     )     ⁢   ℓ     )         ,             Eqn   ⁢           (   78   )                   
and where a chn (f chn (k chn );l) is given by Eqn (40)-Eqn (41) and “⊗” denotes the Kronecker product operation. Similarly, the per-channel and per-subband interference ACM is given by
 
                         R       i   chn     ⁢     i   chn         (     k   chn     )     ≈         M   ADC   2       M   chp       ⁢     f   chp     ⁢   diag   ⁢       {           ❘   &#34;\[LeftBracketingBar]&#34;         H   R     (       f   chn     (     k   chn     )     )       ❘   &#34;\[RightBracketingBar]&#34;       2     ⁢       S       i   R     ⁢     i   R         (       f   chn     (     k   chn     )     )       }         k   sub     ,     =   0           K   sub     -   1           ,           Eqn   ⁢           (   79   )                                   R       i   DoF     ⁢     i   DoF         (     k   sub     )     ≈         M   ADC   2       M   chp       ⁢     f   chp     ⁢         ❘   &#34;\[LeftBracketingBar]&#34;         H   R     (       f   sub     (     k   sub     )     )       ❘   &#34;\[RightBracketingBar]&#34;       2     ⁢     (       P   sub     ⊗       S       i   R     ⁢     i   R         (       f   sub     (     k   sub     )     )       )         ,           Eqn   ⁢           (   80   )                 
where S i     R     i     R   (f) is the M feed ×M feed  matrix power spectral density of i R (t).
 
     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 multi-feed 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. 
     In presence of channel dynamics, the local signal DOA adheres closely to a first-order model over time intervals on the order of 10 seconds in some aspects of the invention, and the individual subbands will experience additional signature blur due to the changing TOA, LOB, and (typically more importantly, due to dependence on changing receiver platform orientation) DOA of the transmitters and receiver. This signature blur is likely to further load the subbands with low-level signature components that will be excised by subsequent linear combining operations. In this regard, the effect of TOA changes is reduced substantively for the multifeed symbol-rate synchronous sub-band channelizer shown in  FIG.  15   , due to the narrower subband channels obtaining as the number of receiver feeds is increased, as shown in  FIG.  16    and  FIG.  17   . 
     Resilient PNT Signal Processing Methods 
     A common set of signal processing methods can be applied to data provided by all of the single-feed and multifeed receiver structures and symbol-rate synchronous channelization and multi-subband formation operations described above. In particular, in-subband linear combining weights approaching the max-SINR solution can be computed using linear-algebraic methods disclosed in the &#39;417 NPA, given any of the following:
         Knowledge of the content of {d T (n sym ;l)} l=1   L     T   , e.g., from navigation data or symbol sequences provided by third-party sources.   Knowledge of a component of {d T (n sym ;l)} l=1   L     T   , e.g., known or periodic pilot signals, as listed in  FIG.  4    for the GPS civil L5, QZSS, Galileo E5AB, E6B E6C, and NAVIC RS BOC signals   Other exploitable structure of {d T (n sym ;l)} l=1   L     T   , e.g., real or BPSK structure, as listed in  FIG.  4    for the GLONASS and GPS L1 legacy signals; QPSK structure, as listed in  FIG.  4    for the Beidou 1B and NAVIC SPS-L5 and SPS-S signals; and constant modulus structure, as listed in  FIG.  4    for all of the GNSS signals.
 
These methods also can be used to estimate the max-SINR obtaining in each subband, thereby providing cross-subband weights for the second combining stage shown in Eqn (64). The multi-subband solution can further be used to estimate the full max-SINR of the received symbol sequences, detect the navigation signals in the environment, determine at least the geo-observables {n TR (l),{tilde over (α)} TR (l)} l=1   L     T    for those signals, and if needed estimate information-bearing components of the transmitted symbol streams {d T (n sym ;l)} l=1   L     T   . Moreover, these techniques can perform these operations without prior knowledge of the signal ranging code or receive frequency signature for the signals, or a search over TOA and FOA is the ranging code is known.
       

     In addition, the &#39;417 NPA discloses means for estimating the fine TOA, FOA Nyquist zone, and DFOA using additional copy-aided parameter estimation algorithms. Lastly, the data model derived above can be used to develop matched-filtering methods, using estimates of a chn (l) given in Eqn (40)-Eqn (41) for single-receiver systems, or defined in Eqn (77) for multi-receiver systems. 
     The signal processing 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 sym  channelized data symbols 
               {       x   DoF     (       n   sym     ;     k   sub       )     }         n   sym     =   0         N   sym     -   1           
are computed for each subband 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 processing algorithms are described in more detail below.
 
Batch Adaptive Processing Procedure
 
     In the batch adaptive processing procedure is implemented by first collecting data over N sym  data symbols, and detecting, extracting, or estimating geo-observables of the MOS-DSSS symbol or navigation sequences directly from that data set. In some aspects, the procedure performs 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&#39;s of the signals (and in particular FOA&#39;s derived from the FOA vectors for those signals) are known, then the procedure can be started at an intermediate point in the processing. 
     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&#39;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). 
       FIG.  18    is a flow diagram of a fully-blind multi-subband signal detection, geo-observable/quality estimation, and symbol estimation batch procedure that can be employed for MOS-DSSS signals with BPSK (or, more generally, real) symbol sequences, e.g., the GLONASS and GPS L1 legacy signal, or LocataLite beacons. Upon receipt of N sym  symbols of data ( 300 ), e.g., N sym  milliseconds of receive data for the GNSS signals depicted in  FIG.  4   , the procedure first channelizes the data and forms it into K sub  subbands using a symbol-rate synchronous channelizer ( 301 ). In some aspects, the data subbands are formed into K sub  data matrices with dimension N sym ×M DoF , e.g., given by 
                         X   DoF     (     k   sub     )       =   Δ       (                 ω   DoF     (   0   )       ⁢       x   DoF   H     (     0   ;     k   sub       )               ⋮                   ω   DoF     (       N   sym     -   1     )       ⁢       x   DoF   H     (         N   sym     -   1     ;     k   sub       )             )       ,           Eqn   ⁢           (   81   )                 
where
 
               {       ω   DoF     (     n   sym     )     }         n   sym     =   0         N   sym     -   1           
is a real data window satisfying
 
     
       
         
           
             
               
                 ∑ 
                 
                   
                     n 
                     sym 
                   
                   = 
                   0 
                 
                 
                   
                     N 
                     sym 
                   
                   - 
                   1 
                 
               
                 
               
                 
                   ω 
                   DoF 
                 
                 ( 
                 
                   n 
                   sym 
                 
                 ) 
               
             
             = 
             1. 
           
         
       
     
     The channelized subband data matrices 
               {       X   DoF     (     k   sub     )     }         k   sub     =   0         K   sub     -   1           
are then whitened ( 302 ), e.g., by performing the QR decomposition (QRD) of each subband matrix, to form whitened subband matrices
 
                 {       Q   DoF     (     k   sub     )     }         k   sub     =   0         K   sub     -   1       .         
The QRD, denoted {Q,R}=QRD(X) for general N×M matrix X, solves
 
                     R   =     chol   ⁢     {       X   H     ⁢   X     }         ,           Eqn   ⁢           (   82   )                                 Q   =       XR     -   1       =     (             q   H     (   0   )             ⋮               q   H     (     N   -   1     )           )         ,           Eqn   ⁢           (   83   )                   
where chol(⋅) is the Cholesky factorization operation, such that X=QR and Q H Q=I M  where I M  M×M identity matrix. The data window is nominally rectangular; however, other windows are recommended if the FOA offset between the navigation signals is small, e.g., for fixed ground beacons or aligned spoofers, or for reception intervals that are long enough to induce substantive differential FOA effects.
 
     For navigation signals in which only a single rail is modulated, e.g., GLONASS or GPS L1 legacy signals, the subbands can be processed using the conjugate self-coherence restoral (CSCORE) method disclosed in the &#39;417 NPA. For the multi-subband method, CSCORE statistics are developed within each subband, and combined to form a cross-subband statistic ( 303 ) that simultaneously detects the signals, and provides an estimate the fine FOA and quality of each signal. In one aspect, the CSCORE algorithm is implemented in each subband for a trial FOA vector 
               α   =       [         α         α     (   1   )             ]     T       ,         
by first computing
 
                   Q   DoF     (     α   ;     k   sub       )     =       Δ   ⁡   (   α   )     ⁢       Q   DoF     (     k   sub     )         ,                   Δ   ⁡   (   α   )       =   Δ       diag   ⁢       {     exp   ⁡   (     j   ⁢   2   ⁢   π   ⁢       ℊ   FOA     (     n   sym     )     ⁢   α     )     }         n   sym     =   0         N   sym     -   1           ,         
where g FOA (n sym )=[n sym  ½n sym   2 ]; forming whitened CSCORE matrix
 
                 R   ^         q   DoF     ⁢     q   DoF   *         (     α   ;     k   sub       )         
given by
 
                               R   ^         q   DoF     ⁢     q   DoF   *         (     α   ;     k   sub       )     =           ∑       n   sym     =   0         N   sym     -   1                 q   DoF     (       n   sym     ;     k   sym       )     ⁢       q   DoF   T     (       n   sym     ;     k   sym       )     ⁢     exp   ⁡   (       -   j     ⁢   4   ⁢   π   ⁢       ℊ   FOA   T     (     n   sym     )     ⁢   α     )                       =           ∑       n   sym     =   0         N   sym     -   1                 q   DoF     (       n   sym     ;   α   ;     k   sub       )     ⁢       q   DoF   T     (       n   sym     ;   α   ;     k   sub       )           ,               =             Q   DoF   H     (     α   ;     k   sub       )     ⁢       Q   DoF   *     (     α   ;     k   sub       )                     Eqn   ⁢           (   84   )                 
for each subband; and determining the dominant mode {λ CSC (α;k sub ),v CSC (α;k sub )} of the singular value decomposition (SVD) of
 
                   R   ^         q   DoF     ⁢     q   DoF   *         (     α   ;     k   sub       )     ,         
e.g., using a power method. The detection statistics
 
               {       λ   CSC     (     α   ;     k   sub       )     }         k   sub     =   0         K   sub     -   1           
are then combined across the subbands to form a full-band CSCORE statistic. In one aspect on the invention, this is accomplished by first computing max-SINR estimate γ CSC (α;k sub )=(1+λ CSC (α;k sub ))/(1−λ CSC (α;k sub )) for each subband, and summing those estimates together for each FOA vector α. In other aspects, nonlinear combining operations, e.g., dictated by maximum-likelihood (ML) estimation arguments are performed. In other aspects, the combining is performed over values of α that are adjusted to account for DTOA between subbands.
 
     This procedure generates a CSCORE spectrum that is parameterized with respect to both FOA (or DTOA) and DFOA. This spectrum is then used to detect the MOS-DSSS signals in the environment, and estimate their FOA (or DTOA) and DFOA geo-observables ( 304 ). Additional joint processing is then performed to further improve quality of the geo-observables and the detection statistic ( 305 ), resulting in optimized DTOA and DFOA estimates {{circumflex over (τ)} CSC   (1) (l),{circumflex over (α)} CSC   (1) (l)} l=1   {circumflex over (L)}     T    for the {circumflex over (L)} T  detected signals, and resulting in optimized despreading weights {v CSC (α(k sub ;l));k sub } l=1   {circumflex over (L)}     T    for each subband, where α(k sub ;l)=[−{circumflex over (τ)} CSC   (1) (l)(f T +f sub (k sub )) {circumflex over (α)} CSC   (1) (l)] T . 
     The optimized despreading weights and estimated FOA vector are then used to despread and despin the symbol sequence for each detected signal, and to demodulate the underlying signal navigation sequence based on further structure of that sequence, e.g., 20:1 NAV bit replication for the GPS L1 legacy signal ( 306 ). As part of this process, the coarse TOA timing is computed to within a NAV bit edge. 
     Given reception of sufficient data, the NAV signal is analyzed to resolve the ±1 ambiguity in the spread signal; detect NAV block boundaries to determine the full coarse TOA; and if needed extract the satellite ephemeres from the NAV signal sequences ( 307 ). In some aspects, this is accomplished by processing of data over multiple reception blocks. In other aspects, this is accomplished by shifting to a continuous update (non-batch) processing mode. 
     In some aspects, the ranging code and (for multifeed receivers) array manifold are downloaded from memory ( 310 ), and used to determine the full TOA, FOA, DOA, and (using estimated and/or on-board orientation data) LOB of the detected signals, and determine a positioning, navigation, and timing (PNT) solution for the receiver ( 308 ). In other aspects, the fine FOA (or TDOA) and DFOA, and transmitter ephemeres are sufficient to determine a positioning, navigation and frequency synchronization solution for the receiver. 
       FIG.  19    is a flow diagram of a multi-subband signal detection and geo-observable/quality estimation, algorithm that can be employed in a resilient PNT data analytics engine (RPNT-AE) aspect of the invention. Upon receipt of an RPNT-AE trigger ( 319 ), e.g., initiated by a primary GNSS processor, or by other processes on board and/or connected to the receiver, the RPNT-AE process is begun ( 320 ). The RPNT-AE processor then tasks a receiver to obtain and store a snapshot of N sym  data symbols, and obtains a time/frequency stamp recording rough (to within receiver clock accuracy) center frequency and start time of the snapshot, and other information available to the receiver, e.g., receiver orientation ( 321 ). The RPNT-AE processor then places the time/frequency stamped data into memory ( 310 ), and sends the stamp to an external resource ( 322   a ) along with a request for third-party data. The RPNT-AE processor then waits until an external resource ( 322   b ) provides the processor with a snapshot of navigation (NAV/CNAV) data covering the time period and frequency band recorded in the time/frequency stamp, and (if needed) a transmitter almanac containing the ephemeres of the transmitters over the time period covering the NAV/CNAV data. Once RPNT-AE processor obtains that data ( 323 ), it places that data in memory ( 310 ) and sends a “Data Ready” flag is sent to the DSP processor tasked with estimating RPNT analytics for the processor. 
     Once the Data Ready flag is obtained from the RPNT-AE processor, the DSP resources obtain the time/frequency stamped received data snapshot and the NAV/CNAV data from memory ( 310 ), and generate RPNT analytics using partially blind methods disclosed in the &#39;417 NPA, e.g., an “FFT-least-squares” algorithm or a single-target or multi-target maximum-likelihood estimator, the DSP processor. Analytics measure here can include the following:
         estimates of observed received parameters of navigation signals, including signal geo-observables usable to obtain PNT solutions for the reception platform, e.g., observed signal frequency-of-arrival (FOA), time-of-arrival (TOA), and observed local direction-of-arrival (DOA) or global line-of-bearing (LOB) in systems employing multifeed receivers;   estimates of received signal quality, e.g., received incident power, and of despread/demodulated navigation sequence quality, e.g., despreader output signal-to-interference-and-noise ratio (SINR); and   measures of accuracy of parameter and signal quality estimates.
 
RPNT analytic measurement methods include “partially-blind” methods that do not require knowledge of the ranging code for the navigation signals, or spatial/polarization array-manifold data (measurement of cross-feed spatial/polarization signatures as a function of DOA) for the receiver, and “copy-aided” methods that may require one or both of the ranging codes or the array manifold to provide more complete RPNT analytics. In the aspects disclosed herein, RPNT analytics are further computed over multiple subbands.
       

     Once the RPNT analytics have been obtained, the RPNT-AE process is ended ( 325 ), and the RPNT analytics are reported to the resources requesting the analytics. 
       FIG.  20    is a block diagram of a multi-subband pilot-exploiting RPNT (PE-RPNT) processor that can detect and estimate geo-observables and quality of MOS-DSSS signals in which known or repeating pilots transmitted as part of the MOS-DSSS system sequence. The method reported herein can also be used to implement the operations performed by the RPNT-AE DSP processor ( 324 ), using the full MOS-DSSS symbol sequence provided in that system. 
     Upon receipt of N sym  symbols of data ( 300 ), e.g., N sym  milliseconds of receive data for the GNSS signals depicted in  FIG.  4   , the procedure first channelizes the data and forms it into K sub  subbands using a symbol-rate synchronous channelizer ( 301 ). In some aspects, the data subbands are formed into K sub  data matrices with dimension N sym ×M DoF , e.g., given by Eqn (81). The channelized subband data matrices 
               {       X   DoF     (     k   sub     )     }         k   sub     =   0         K   sub     -   1           
are then whitened ( 302 ), e.g., by performing the QR decomposition (QRD) of each subband matrix, to form whitened subband matrices
 
     
       
         
           
             
               
                 { 
                 
                   
                     Q 
                     DoF 
                   
                   ( 
                   
                     k 
                     sub 
                   
                   ) 
                 
                 } 
               
               
                 
                   k 
                   sub 
                 
                 = 
                 0 
               
               
                 
                   K 
                   sub 
                 
                 - 
                 1 
               
             
             . 
           
         
       
     
     The processor then retrieves a single period of repeating pilot data 
               {     p   ⁡   (     n   sym     )     }         n   sym     =   0         N   pilot     -   1           
from memory ( 310 ), where N pilot  is the period of the pilot sequence; computes pilot statistics over a set of trial FOA vectors
 
             α   =       [         α         α     (   1   )             ]     T           
trial pilot delay m pilot ∈{0, . . . , N pilot −1} in each subband, given by least-squares (LS) whitened linear combiner weights
 
                             u   LS     (     α   ,       m   pilot     ;     k   sub         )     =           ∑       n   sym     =   0       N   sym                 q   DoF     (       n   sym     ;     k   sub       )     ⁢       d   LS   *     (         n   sym     ;   α     ,     m   pilot       )                       =             Q   DoF   H     (     k   sub     )     ⁢       d   LS     (     α   ,     m   pilot       )         ,                 Eqn   ⁢           (   85   )                           d   LS     (         n   sym     ;   α     ,     m   sym       )     =           ω   DoF     (     n   sym     )       ⁢     p   ⁡   (       (       n   sym     -     m   pilot       )     ⁢   mod   ⁢     N   pilot       )     ⁢     e     j   ⁢   2   ⁢   π   ⁢       ℊ   FOA     (     n   sym     )     ⁢   α                             	           d   LS     (     α   ,     m   sym       )     =     (             d   LS   *     (       0   ;   α     ,     m   sym       )             ⋮               d   LS   *     (           N   sym     -   1     ;   α     ,     m   sym       )           )       ,             Eqn   ⁢           (   86   )                   
where g FOA (n sym )=[n sym  ½n sym   2 ], and in-band LS SINR estimate,
 
                         γ   LS     (     α   ,       m   pilot     ;     k   sub         )     =                u   LS     (     α   ,       m   pilot     ;     k   sub         )          2   2       1   -              u   LS     (     α   ,       m   pilot     ;     k   sub         )          2   2           ;           Eqn   ⁢           (   87   )                 
and combines the in-subband least-square SINR estimates to created a multi-subband quality statistic as a function of FOA vector and pilot delay ( 333 ). In one aspect, the in-subband LS SINR estimates are combined in accordance with a maximum-likelihood estimator, yielding multi-subband quality statistic
 
                         S   ML     (     α   ,     m   pilot       )     =       ∑       k   sub     =   0         K   sub     -   1             log   ⁡   (     1   +       γ   LS     (     α   ,       m   pilot     ;     k   sub         )       )         ,           Eqn   ⁢           (   88   )                 
which is maximized at modulo-N pilot  coarse TOA&#39;s and fine FOA vectors possessed by the MOS-DSSS signals in the receivers field of view. In another aspect, the in-subband LS SINR estimates are combined in accordance with Eqn (65), yielding multi-subband ML quality statistic
 
                         S   LS     (     α   ,     m   pilot       )     =       ∑       k   sub     =   0         K   sub     -   1               γ   LS     (     α   ,       m   pilot     ;     k   sub         )         ,           Eqn   ⁢           (   89   )                 
which is also maximized at modulo-N pilot  coarse TOA&#39;s and fine FOA vectors possessed by the MOS-DSSS signals in the receiver&#39;s field of view. In other aspects, the effect of DTOA across frequency subbands is taken into account, e.g., resulting in DTOA-corrected multi-subband LS quality statistic
 
                         S   LS     (       τ     (     1   ,   2     )       ,     m   pilot       )     =       ∑       k   sub     =   0         K   sub     -   1               γ   LS     (       α   ⁡   (     k   sub     )     ,       m   pilot     ;     k   sub         )         ,           Eqn   ⁢           (   90   )                                 τ     (     1   ,   2     )       =     (           τ     (   1   )                   τ   ~       (   2   )             )       ,           Eqn   ⁢           (   91   )                                 α   ⁡   (     k   sub     )     =     -     (       f   T     +       f   sub     (     k   sub     )       )     ⁢     T   sym     ⁢     τ     (     1   ,   2     )           ,             Eqn   ⁢           (   92   )                   
where {tilde over (τ)} (2) =τ (2) T sym  is the symbol-normalized differential DTOA. In addition, for rectangular time windows, the in-band LS SINR estimate can be converted to unbiased max-SINR estimate
 
                           γ   ^       max   -   SINR       (     α   ,       m   pilot     ;     k   sub         )     =         (     1   -         M   DoF     +   1       N   sym         )     ⁢       γ   LS     (     α   ,       m   pilot     ;     k   sub         )       -     M   DoF         ,           Eqn   ⁢           (   93   )                 
which can greatly improve visible quality of S LS (α,m pilot ) (or S LS (τ (1,2) ;m pilot )) without affecting performance of the statistic.
 
     A number of means can be used to minimize efficiency of the in-subband LS SINR&#39;s and multi-subband quality statistics. In particular, FFT-based methods can be used to efficiently compute γ LS (α,m pilot ;k sub ) or γ LS (α(k sub )m pilot ;k sub ) with fine accuracy. At this stage of processing, the whitened LS weights given in Eqn (85) need not be computed, further reducing complexity and memory requirements of the overall processor. 
     The multi-subband quality statistic is then analyzed to detect the pilot-bearing MOS-DSSS signals in the environment, determine their coarse TOA to modulo-N pilot  accuracy, and determine their fine FOA or DTOA vector ( 334 ), resulting in {circumflex over (L)} T  peak detections and FOA-TOA locations {α(l),m pilot (l)} l=1   {circumflex over (L)}     T    or DTOA-TOA locations {τ (1,2) (l),m pilot (l)} l=1   {circumflex over (L)}     T    of those peaks. Iterative refinement of the DTOA-TOA locations {τ (1,2) (l),m pilot (l)} l=1   {circumflex over (L)}     T    (where the FOA vector is converted to a DTOA vector) is then performed ( 335 ), e.g., using Newton search methods disclosed in the &#39;417 NPA, to determine the fine DTOA vector (from which FOA can be derived) and optimize the multi-subband quality statistic. 
     Using the optimized DTOA vectors and modulo-N pilot  coarse TOA&#39;s {τ (1,2) (l),m pilot (l)} l=1   {circumflex over (L)}     T   , LS combiner weights {γ LS (k sub ;l),u LS (k sub ;l)} are then computed using formula 
                         u   LS     (       k   sub     ;   ℓ     )     =       ∑       n   sym     =   0       N   sym                 q   DoF     (       n   sym     ;     k   sub       )     ⁢       d   LS   *     (         n   sym     ;     α   ⁡   (       k   sub     ;   ℓ     )       ,       m   pilot     (   ℓ   )       )           ,           Eqn   ⁢           (   94   )                             	         α   ⁡   (       k   sub     ;   ℓ     )     =     -     (       f   T     +       f   sub     (     k   sub     )       )     ⁢     T   sym     ⁢       τ     (     1   ,   2     )       (   ℓ   )         ,               Eqn   ⁢           (   95   )                               	           γ   LS     (       k   sub     ;   ℓ     )     =                u   LS     (       k   sub     ;   ℓ     )          2   2       1   -              u   LS     (       k   sub     ;   ℓ     )          2   2           ,             Eqn   ⁢           (   96   )                 
and cross-subband combiner gains are computed based on Eqn (64), i.e., by setting g LS (k sub ;l)=1+γ LS (k sub ;l) The full multi-subband despread, despun and partially time-synchronized symbol sequence is then estimated for each peak using formula
 
                   	             d   ^     T     (       n   sym     ;   ℓ     )     =       ∑       k   sub     =   0         K   sub     -   1                 ℊ   LS   *     (       k   sub     ;   ℓ     )     ⁢         d   ^     T     (       n   sym     ;     k   sub     ;   ℓ     )           ,             Eqn   ⁢           (   97   )                                   d   ^     T     (       n   sym     ;     k   sub     ;   ℓ     )     =       (         u   LS   H     (     k   sub     )     ⁢       q   DoF     (         n   sym     +       m   pilot     (   ℓ   )       ;   ℓ     )       )     ⁢       e     -   j   ⁢   2   ⁢   π   ⁢       ℊ   FOA     (     n   sym     )     ⁢     α   ⁡   (       k   sub     ;   ℓ     )         .               Eqn   ⁢           (   98   )                 
The navigation-bearing component of the symbol sequence (CNAV for the GPS L5 signal) is then demodulated using information about the structure of the navigation data ( 346 ).
 
     Given reception of sufficient data, the CNAV signal is analyzed to determine the full coarse TOA, and if needed extract the satellite ephemeres from the NAV signal sequences ( 347 ). In some aspects, this is accomplished by processing of data over multiple reception blocks. In other aspects, this is accomplished by shifting to a continuous update (non-batch) processing mode. 
     In some aspects, the ranging code and (for multifeed receivers) array manifold are downloaded from memory ( 310 ), and used to determine the full TOA, FOA, DOA, and (using estimated and/or on-board orientation data) LOB of the detected signals, and determine a positioning, navigation, and timing (PNT) solution for the receiver ( 308 ). In other aspects, the fine FOA (or TDOA) and DFOA, and transmitter ephemeres are sufficient to determine a positioning, navigation and frequency synchronization solution for the receiver. 
       FIG.  21    is a flow diagram of a local-maxima search procedure according to an aspect of the disclosure, which is applicable to searches over FOA and DFOA, over reception intervals in which DFOA can have a substantive effect on detection sensitivity. The approach assumes the use of an FFT-based coarse search procedure, in which α is defined over search grid 
                         α   ~     (     k   bin     )     =       k         bin       /     K   bin         ,       k   bin     =   0     ,   …         ,       K   bin     -   1     ,           Eqn   ⁢           (   99   )                                           α   ~     TR     (   1   )       (     k   tile     )     =       α   ~     TR     (   1   )           ❘   &#34;\[RightBracketingBar]&#34;       max     ⁢     (           2   ⁢     k   tile       +   1       K   tile       -   1     )       ,       k   tile     =   0     ,   …         ,       K   tile     -   1     ,             Eqn   ⁢           (   100   )                   
where K bin  is the number of FOA bins searched over, and K tile  is the number of DFOA tiles searched over. Then the multi-subband quality estimate can be given by 3-dimensional surface
 
                         S   LS     (       k   bin     ,     k   tile     ,     m   pilot       )     =       ∑       k   sub     =   0         K   sub     -   1               γ   LS     (       k   bin     ,     k   tile     ,       m   pilot     ;     k   sub         )         ,           Eqn   ⁢           (   101   )                                   γ   LS     (       k   bin     ,     k   tile     ,       m   pilot     ;     k   sub         )     =                u   LS     (       k   bin     ,     k   tile     ,       m   pilot     ;     k   sub         )          2   2       1   -              u   LS     (       k   bin     ,     k   tile     ,       m   pilot     ;     k   sub         )          2   2           ,           Eqn   ⁢           (   102   )                                       u   LS     (       k   bin     ,     k   tile     ,       m   pilot     ;     k   sub         )     =           ∑       n   sym     =   0         N   sym     -   1                 q   DoF     (       n   sym     ;     k   sub       )     ⁢       d   LS   *     (         n   sym     ;     k   tile       ,     m   pilot       )     ⁢     e     -   j   ⁢   2   ⁢   π   ⁢     k   bin     ⁢     n   sym     /     K   bin                           =         DFT   ⁢     {         q   DoF     (       n   sym     ;     k   sub       )     ⁢       d   LS   *     (         n   sym     ;     k   tile       ,     m   pilot       )       }         ,                 Eqn   ⁢           (   103   )                                   d   LS     (         n   sym     ;     k   tile       ,     m   pilot       )     =       e     j   ⁢   π   ⁢         α   ~       (   1   )       (     k   tile     )     ⁢     n   sym   2         ⁢     p   ⁡   (       (       n   sym     -     m   pilot       )     ⁢   mod   ⁢     N   pilot       )         ,           Eqn   ⁢           (   104   )                 
which can be implemented using K tile N pilot  DFT operations.
 
     Returning to  FIG.  21   , once the search is started ( 350 ), the number of DFOA search tiles K tile  first is determined based on the expected range in DFOA and the reception time. Then an initial DFOA search tile is chosen ( 351 ), and a DFT-based search is conducted over that tile using Eqn (101)-Eqn (104), to compute a multi-subband search spectrum for that tile, and determine peaks in that tile ( 352 ). 
     Once all of the DFOA tiles have been searched, and peaks have been determined, the detected peaks are thinned to a significant number of peaks based on peak value ( 353 ). Then the peak values are iterivaley optimized using local optimization methods, e.g., quadratic fit to the peak maximum, followed by Newton or Gauss-Newton search methods ( 354 ). In this case, m pilot  is defined over an integer set of values, hence it is not iteratively optimized using this procedure. 
     In one aspect of the invention, the local optimization is performed to optimize 
                         S   LS     (   ℓ   )     =       ∑       k   sub     =   0         K   sub     -   1               γ   LS     (       k   sub     ;   ℓ     )         ,           Eqn   ⁢           (   105   )                 
where γ LS (k sub ;l) is given by Eqn (94)-Eqn (96), and where {τ (1,2) (l),m pilot (l)} l=1   {circumflex over (L)}     T    is the parameter set being optimized. In this case, τ (1,2) (l) is initialized using the optimal value of α(l) determined during the FFT-based search procedure.
 
     Once optimal values have been determined, ancillary RPNT statistics are computed at those optimal search locations ( 355 ). The search is then ended ( 356 ). 
       FIGS.  22 A,  22 B, and  22 C  and  FIGS.  23 A,  23 B, and  23 C  depict DFOA search sensitivity as a function of the number of DFOA search tiles, for a ±3 Hz/s DFOA range, consistent with DFOA&#39;s expected for a MEO satellite constellation.  FIGS.  22 A,  22 B, and  22 C  show search sensitivity relative to 0 Hz/s for a rectangular time window, and for a ½ second ( FIG.  22 A ), 1 second ( FIG.  22 B ), and 2 second ( FIG.  22 C ) reception interval. As shown in  FIGS.  22 A,  22 B, and  22 C , the number of search tiles needed to minimize sensitivity loss grows rapidly with reception time.  FIGS.  23 A,  23 B, and  23 C  show equivalent search sensitivity for a Nuttall time window, which has a time-bandwidth product (TBP) that is roughly half the value of the rectangular window. The Nuttall window can reduce the number of tiles required to perform an effective search at a 2 second reception time. However, this is offset by the lower TBP of the window. 
     While this invention is susceptible of instantiation in many different forms, there are shown in the drawings and described in detail in the text of Provisional Appl. No. 62/773,589, incorporated herein by reference, and Provisional Appl. No. 62/773,605, incorporated herein by reference, several specific aspects of the invention, with the understanding that the present disclosure is to be considered as an exemplification of the principles of the invention and is not intended to limit the invention to the aspects illustrated. 
     In an aspects applicable to all of the approaches above, the invention obtains snapshots of baseband navigation data covering the time interval of data collected by the receiver, and symbol-synchronously channelized by the invention, and uses that baseband navigation data to implement non-fully-blind demodulation algorithms. The resultant aspects can provide extreme high precision of FOA, TOA, and DFOA/DTOA drift estimates; assess integrity of data collected by the host platform; and provide other functions of interest to the user. When coupled with a communication channel allowing the receive data to be transported to a central processor, the invention can also allow implementation of all functions at off-line resources, thereby eliminating all DSP complexity associated with the algorithms. The aspects can also be used to implement signal cancellation algorithms that detect signals under the known navigation signals, e.g., for purposes of spoofer and jammer detection. 
     Any reception operation used in the invention can be implemented using any of the set of one or more dedicated receivers and software defined radios (SDR) either separate from or integrated with antennas, amplifiers, mixers, filters, analog-to-digital converters (ADC&#39;s) and signal processing gear. 
     Operations Processing used in each of the inventions above can be implemented in any combination of hardware and software, from special-purpose hardware including any of application-specific integrated circuits (ASIC&#39;s) and field-programmable gate arrays (FPGA&#39;s); firmware instructions in a lesser-specialized set of hardware; embedded digital signal processors (e.g. Texas Instrument or Advanced Risc Machine DSP&#39;s); graphical processing units (GPU&#39;s); vector, polynomic, quantum, and other processors; and in any combination or sole use of serial or parallel processing; and on general-purpose computers using software instructions. 
     Operations Processing used in each of the inventions above can be further implemented using any set of resources that are on-board, locally accessible to, and remotely accessible by the receiver after transport of the data and instructions to be processed by any of a single computer, server, and set of servers, and then directed onwards, using any number of wired or wireless means for such transport. 
     Some of the above-described functions may be composed of instructions, or depend upon and use data, that are stored on storage media (e.g., computer-readable medium). Some of the above-described functions may be comprised in EEPROMs, ASICs, or other combinations of digital circuitry for digital signal processing, connecting and operating with the adaptive processor. The instructions and/or data may be retrieved and executed by the adaptive processor. Some examples of storage media are memory devices, tapes, disks, and the like. The instructions are operational when executed by the adaptive processor to direct the adaptive processor to operate in accord with the invention; and the data is used when it forms part of any instruction or result therefrom. 
     The terms “computer-readable storage medium” and “computer-readable storage media” as used herein refer to any medium or media that participate in providing instructions to a CPU for execution. Such media can take many forms, including, but not limited to, non-volatile (also known as ‘static’ or ‘long-term’) media, volatile media and transmission media. Non-volatile media include, for example, one or more optical or magnetic disks, such as a fixed disk, or a hard drive. Volatile media include dynamic memory, such as system RAM or transmission or bus ‘buffers’. Common forms of computer-readable media include, for example, a floppy disk, a flexible disk, a hard disk, magnetic tape, any other magnetic medium, a CD-ROM disk, digital video disk (DVD), any other optical medium, any other physical medium with patterns of marks or holes. 
     Memory, as used herein when referencing to computers, is the functional hardware that for the period of use retains a specific structure which can be and is used by the computer to represent the coding, whether data or instruction, which the computer uses to perform its function. Memory thus can be volatile or static, and be any of a RAM, a PROM, an EPROM, an EEPROM, a FLASH-EPROM, any other memory chip or cartridge, a carrier wave, or any other medium from which a computer can read data, instructions, or both. 
     I/O, or ‘input/output’, is any means whereby the computer can exchange information with the world external to the computer. This can include a wired, wireless, acoustic, infrared, or other communications link (including specifically voice or data telephony); a keyboard, tablet, camera, video input, audio input, pen, or other sensor; and a display (2D or 3D, plasma, LED, CRT, tactile, or audio). That which allows another device, or a human, to interact with and exchange data with, or control and command, a computer, is an I/O device, without which any computer (or human) is essentially in a solipsistic state. 
     While this invention has been described in reference to illustrative aspects of the invention, this description is not to be construed in a limiting sense. Various modifications and combinations of the illustrative aspects as well as other aspects of the invention will be apparent to those skilled in the art upon referencing this disclosure. It is therefore intended this disclosure encompass any such modifications or aspects. 
     The scope of this invention includes any combination of the elements from different aspects disclosed in this specification, and is not limited to the specifics of any of the aspects mentioned above. Individual user configurations and aspects of this invention may contain all, or less than all, of the elements disclosed in the specification according to the needs and desires of that user. The claims stated herein should be read as including those elements which are not necessary to the invention yet are in the prior art and are necessary to the overall function of that particular claim, and should be read as including, to the maximum extent permissible by law, known functional equivalents to the elements disclosed in the specification, even though those functional equivalents are not exhaustively detailed herein. 
     Although the present invention has been described chiefly in terms of the specific aspects of the invention, it is to be understood that the disclosure is not to be interpreted as limiting. Various alterations and modifications will no doubt become apparent to those skilled in the art after having read the above disclosure. Such modifications may involve other features which are already known in the design, manufacture and use of wireless electromagnetic communications networks, systems and MIMO networks and systems therefore, and which may be used instead of or in addition to features already described herein. The algorithms and equations herein are not limiting but instructive of aspects of the invention, and variations which are readily derived through programming or mathematical transformations which are standard or known to the appropriate art are not excluded by omission. Accordingly, it is intended that the appended claims are interpreted as covering all alterations and modifications as fall within the true spirit and scope of the invention in light of the prior art. 
     Additionally, although claims have been formulated in this application to particular combinations of elements, it should be understood that the scope of the disclosure of the present application also includes any single novel element or any novel combination of elements disclosed herein, either explicitly or implicitly, whether or not it relates to the same invention as presently claimed in any claim and whether or not it mitigates any or all of the same technical problems as does the present invention. The applicants hereby give notice that new claims may be formulated to such features and/or combinations of such features during the prosecution of the present application or of any further application derived therefrom.