Add acquisition support for Glonass

Depends on
#26

Author: Pasi Miettinen

peregrine/acquisition.py

@@ -20,6 +20,8 @@
 import defaults
 
 from include.generateCAcode import caCodes
+from include.glo_ca_code import value as glo_ca_code
+from peregrine.gps_constants import L1CA, GLO_L1, glo_l1_step
 
 import logging
 logger = logging.getLogger(__name__)
@@ -188,13 +190,15 @@ def interpolate(self, S_0, S_1, S_2, interpolation='gaussian'):
 
     **Parabolic interpolation:**
 
-    .. math:: \Delta = \\frac{1}{2} \\frac{S[k+1] - S[k-1]}{2S[k] - S[k-1] - S[k+1]}
+    .. math:: \Delta = \\frac{1}{2} \\frac{S[k+1] -
+                         S[k-1]}{2S[k] - S[k-1] - S[k+1]}
 
     Where :math:`S[n]` is the magnitude of FFT bin :math:`n`.
 
     **Gaussian interpolation:**
 
-    .. math:: \Delta = \\frac{1}{2} \\frac{\ln(S[k+1]) - \ln(S[k-1])}{2\ln(S[k]) - \ln(S[k-1]) - \ln(S[k+1])}
+    .. math:: \Delta = \\frac{1}{2} \\frac{\ln(S[k+1]) -
+                         \ln(S[k-1])}{2\ln(S[k]) - \ln(S[k-1]) - \ln(S[k+1])}
 
     The Gaussian interpolation method gives better results, especially when
     used with a Gaussian window function, at the expense of computational
@@ -363,7 +367,9 @@ def find_peak(self, freqs, results, interpolation='gaussian'):
     return (code_phase, freq, snr)
 
   def acquisition(self,
-                  prns=range(32),
+                  bandcode=L1CA,
+                  prns=xrange(32),
+                  channels=[x - 7 for x in xrange(14)],
                   doppler_priors=None,
                   doppler_search=7000,
                   doppler_step=None,
@@ -372,10 +378,10 @@ def acquisition(self,
                   multi=True
                   ):
     """
-    Perform an acquisition for a given list of PRNs.
+    Perform an acquisition for a given list of PRNs/channels.
 
-    Perform an acquisition for a given list of PRNs across a range of Doppler
-    frequencies.
+    Perform an acquisition for a given list of PRNs/channels across a range of
+    Doppler frequencies.
 
     This function returns :class:`AcquisitionResult` objects containing the
     location of the acquisition peak for PRNs that have an acquisition
@@ -387,8 +393,13 @@ def acquisition(self,
 
     Parameters
     ----------
+    bandcode : optional
+      String defining the acquisition code. Default: L1CA
+      choices: L1CA, GLO_L1 (in gps_constants.py)
     prns : iterable, optional
       List of PRNs to acquire. Default: 0..31 (0-indexed)
+    channels : iterable, optional
+      List of channels to acquire. Default: -7..6
     doppler_prior: list of floats, optional
       List of expected Doppler frequencies in Hz (one per PRN).  Search will be
       centered about these.  If None, will search around 0 for all PRNs.
@@ -406,10 +417,15 @@ def acquisition(self,
     Returns
     -------
     out : [AcquisitionResult]
-      A list of :class:`AcquisitionResult` objects, one per PRN in `prns`.
+      A list of :class:`AcquisitionResult` objects, one per PRN in `prns` or
+      channel in 'channels'.
 
     """
-    logger.info("Acquisition starting")
+    if bandcode != L1CA and bandcode != GLO_L1:
+      logger.critical("Unkown/Unsupported code " + bandcode)
+      return
+
+    logger.info("Acquisition starting for " + bandcode)
     from peregrine.parallel_processing import parmap
 
     # If the Doppler step is not specified, compute it from the coarse
@@ -421,9 +437,6 @@ def acquisition(self,
       # magnitude.
       doppler_step = self.sampling_freq / self.n_integrate
 
-    if doppler_priors is None:
-      doppler_priors = np.zeros_like(prns)
-
     if progress_bar_output == 'stdout':
       show_progress = True
       progress_fd = sys.stdout
@@ -439,33 +452,55 @@ def acquisition(self,
       show_progress = False
       logger.warning("show_progress = True but progressbar module not found.")
 
+    if bandcode == L1CA:
+      input_len = len(prns)
+      offset = 1
+      pb_attr = progressbar.Attribute('prn', '(PRN: %02d)', '(PRN --)')
+      if doppler_priors is None:
+        doppler_priors = np.zeros_like(prns)
+    else:
+      input_len = len(channels)
+      offset = 0
+      pb_attr = progressbar.Attribute('ch', '(CH: %02d)', '(CH --)')
+      if doppler_priors is None:
+        doppler_priors = np.zeros_like(channels)
+
     # Setup our progress bar if we need it
     if show_progress and not multi:
       widgets = ['  Acquisition ',
-                 progressbar.Attribute('prn', '(PRN: %02d)', '(PRN --)'), ' ',
+                 pb_attr, ' ',
                  progressbar.Percentage(), ' ',
                  progressbar.ETA(), ' ',
                  progressbar.Bar()]
       pbar = progressbar.ProgressBar(widgets=widgets,
                                      maxval=int(len(prns) *
-                                                (2 * doppler_search / doppler_step + 1)),
+                                       (2 * doppler_search / doppler_step + 1)),
                                      fd=progress_fd)
       pbar.start()
     else:
       pbar = None
 
     def do_acq(n):
-      prn = prns[n]
+      if bandcode == L1CA:
+        obj = prns[n]
+        code = caCodes[obj]
+        int_f = self.IF
+        attr = {'prn': obj + 1}
+      else:
+        obj = channels[n]
+        code = glo_ca_code
+        int_f = self.IF + obj * glo_l1_step
+        attr = {'ch': obj}
       doppler_prior = doppler_priors[n]
       freqs = np.arange(doppler_prior - doppler_search,
-                        doppler_prior + doppler_search, doppler_step) + self.IF
+                        doppler_prior + doppler_search, doppler_step) + int_f
       if pbar:
         def progress_callback(freq_num, num_freqs):
-          pbar.update(n * len(freqs) + freq_num, attr={'prn': prn + 1})
+          pbar.update(n * len(freqs) + freq_num, attr=attr)
       else:
         progress_callback = None
 
-      coarse_results = self.acquire(caCodes[prn], freqs,
+      coarse_results = self.acquire(code, freqs,
                                     progress_callback=progress_callback)
 
       code_phase, carr_freq, snr = self.find_peak(freqs, coarse_results,
@@ -478,13 +513,23 @@ def progress_callback(freq_num, num_freqs):
         status = 'A'
 
       # Save properties of the detected satellite signal
-      acq_result = AcquisitionResult(prn,
-                                     carr_freq,
-                                     carr_freq - self.IF,
-                                     code_phase,
-                                     snr,
-                                     status,
-                                     'l1ca')
+      if bandcode == L1CA:
+        acq_result = AcquisitionResult(obj,
+                                       carr_freq,
+                                       carr_freq - self.IF,
+                                       code_phase,
+                                       snr,
+                                       status,
+                                       L1CA)
+      else:
+        acq_result = GloAcquisitionResult(obj,
+                                          carr_freq,
+                                          carr_freq -
+                                            self.IF - obj * glo_l1_step,
+                                          code_phase,
+                                          snr,
+                                          status,
+                                          GLO_L1)
 
       # If the acquisition was successful, log it
       if (snr > threshold):
@@ -494,9 +539,9 @@ def progress_callback(freq_num, num_freqs):
 
     if multi:
       acq_results = parmap(
-          do_acq, range(len(prns)), show_progress=show_progress)
+        do_acq, xrange(input_len), show_progress=show_progress)
     else:
-      acq_results = map(do_acq, range(len(prns)))
+      acq_results = map(do_acq, xrange(input_len))
 
     # Acquisition is finished
 
@@ -505,9 +550,11 @@ def progress_callback(freq_num, num_freqs):
       pbar.finish()
 
     logger.info("Acquisition finished")
-    acquired_prns = [ar.prn + 1 for ar in acq_results if ar.status == 'A']
-    logger.info("Acquired %d satellites, PRNs: %s.",
-                len(acquired_prns), acquired_prns)
+    acq = [ar.prn + offset for ar in acq_results if ar.status == 'A']
+    if bandcode == L1CA:
+      logger.info("Acquired %d satellites, PRNs: %s.", len(acq), acq)
+    else:
+      logger.info("Acquired %d channels: %s.", len(acq), acq)
 
     return acq_results
 
@@ -524,7 +571,7 @@ def save_wisdom(self, wisdom_file=DEFAULT_WISDOM_FILE):
           pyfftw.export_wisdom(), f, protocol=cPickle.HIGHEST_PROTOCOL)
 
 
-class AcquisitionResult:
+class AcquisitionResult(object):
   """
   Stores the acquisition parameters of a single satellite.
 
@@ -553,7 +600,7 @@ class AcquisitionResult:
   """
 
   __slots__ = ('prn', 'carr_freq', 'doppler',
-               'code_phase', 'snr', 'status', 'signal')
+               'code_phase', 'snr', 'status', 'signal', 'sample_index')
 
   def __init__(self, prn, carr_freq, doppler, code_phase, snr, status, signal,
                sample_index=0):
@@ -608,6 +655,20 @@ def _equal(self, other):
     return True
 
 
+class GloAcquisitionResult(AcquisitionResult):
+
+  def __init__(self, channel, carr_freq, doppler, code_phase, snr, status,
+               signal, sample_index=0):
+    super(GloAcquisitionResult, self).__init__(channel, carr_freq, doppler,
+                                               code_phase, snr, status,
+                                               signal, sample_index)
+
+  def __str__(self):
+    return "CH %2d (%s) SNR %6.2f @ CP %6.1f, %+8.2f Hz %s" % \
+        (self.prn, self.signal, self.snr, self.code_phase, self.doppler,
+          self.status)
+
+
 def save_acq_results(filename, acq_results):
   """
   Save a set of acquisition results to a file.
@@ -669,4 +730,5 @@ def print_scores(acq_results, pred, pred_dopp=None):
 
   print "Found %d of %d, mean doppler error = %+5.0f Hz, mean abs err = %4.0f Hz, worst = %+5.0f Hz"\
         % (n_match, len(pred),
-           sum_dopp_err / max(1, n_match), sum_abs_dopp_err / max(1, n_match), worst_dopp_err)
+           sum_dopp_err / max(1, n_match), sum_abs_dopp_err /
+           max(1, n_match), worst_dopp_err)

peregrine/defaults.py

@@ -11,13 +11,18 @@
 ms_to_track = 37 * 1e3
 skip_samples = 1000
 file_format = 'piksi'
-processing_block_size = 20e6 # [samples]
+processing_block_size = 20e6  # [samples]
 
 chipping_rate = 1.023e6  # Hz
 code_length = 1023  # chips
 
 code_period = code_length / chipping_rate
 
+glo_chipping_rate = 0.511e6  # Hz
+glo_code_length = 511  # chips
+
+glo_code_period = glo_code_length / glo_chipping_rate
+
 # original
 sample_channel_GPS_L1 = 0
 sample_channel_GPS_L2 = 1
@@ -80,12 +85,14 @@
 freq_profile_normal_rate = {
     'GPS_L1_IF': 14.58e6,
     'GPS_L2_IF': 7.4e6,
+    'GLO_L1_IF': 12e6,
     'sampling_freq': 24.84375e6}
 
-# 'normal_rate' frequencies profile
+# 'high_rate' frequencies profile
 freq_profile_high_rate = {
     'GPS_L1_IF': freq_profile_normal_rate['GPS_L1_IF'],
     'GPS_L2_IF': freq_profile_normal_rate['GPS_L2_IF'],
+    'GLO_L1_IF': freq_profile_normal_rate['GLO_L1_IF'],
     'sampling_freq': 99.375e6}
 
 L1CA_CHANNEL_BANDWIDTH_HZ = 1000

peregrine/gps_constants.py

@@ -2,23 +2,33 @@
 
 # Some fundamental constants have specific numeric definitions to ensure
 # consistent results in curve fits:
-c = 2.99792458e8 # m/s
+c = 2.99792458e8  # m/s
 pi = 3.1415926535898
 
 # Physical parameters of the Earth
-earth_gm = 3.986005e14 # m^3/s^2 (WGS84 earth's gravitational constant)
-omegae_dot = 7.2921151467e-005 # rad/s (WGS84 earth rotation rate)
+earth_gm = 3.986005e14  # m^3/s^2 (WGS84 earth's gravitational constant)
+omegae_dot = 7.2921151467e-005  # rad/s (WGS84 earth rotation rate)
 
 # GPS system parameters:
-l1 = 1.57542e9 # Hz
-l2 = 1.22760e9 # Hz
+l1 = 1.57542e9  # Hz
+l2 = 1.22760e9  # Hz
 chips_per_code = 1023
-chip_rate = 1.023e6 # Hz
-nominal_range = 26000e3 # m
+chip_rate = 1.023e6  # Hz
+nominal_range = 26000e3  # m
+
+# GLO system parameters
+glo_l1 = 1.602e9  # Hz
+glo_l2 = 1.246e9  # Hz
+glo_chips_per_code = 511
+glo_chip_rate = 0.511e6  # Hz
+glo_l1_step = 0.5625e6  # Hz
 
 # Useful derived quantities:
 code_period = chips_per_code / chip_rate
 code_wavelength = code_period * c
+glo_code_period = glo_chips_per_code / glo_chip_rate
+glo_code_wavelength = glo_code_period * c
 
 L1CA = 'l1ca'
 L2C = 'l2c'
+GLO_L1 = 'glo_l1'

peregrine/initSettings.py

@@ -15,9 +15,12 @@ def __init__(self, freq_profile):
     self.skipNumberOfBytes    = 0                       # Skip bytes in sample file before loading samples for acquisition (bytes)
     self.L1_IF                = freq_profile['GPS_L1_IF']   # L1 intermediate frequency of signal in sample file (Hz)
     self.L2_IF                = freq_profile['GPS_L2_IF']   # L2 intermediate frequency of signal in sample file (Hz)
+    self.GLO_L1_IF            = freq_profile['GLO_L1_IF']   # GLO L1 intermediate frequency of signal in sample file (Hz)
     self.samplingFreq         = freq_profile['sampling_freq'] # Sampling frequency of sample file (Hz)
     self.codeFreqBasis        = defaults.chipping_rate  # Frequency of chipping code (Hz)
     self.codeLength           = defaults.code_length    # Length of chipping code (chips)
+    self.gloCodeFreqBasis     = defaults.glo_chipping_rate  # GLO frequency of chipping code (Hz)
+    self.gloCodeLength        = defaults.glo_code_length    # GLO length of chipping code (chips)
     self.acqThreshold         = 21.0                    # SNR (unitless)
     self.acqSanityCheck       = True                    # Check for sats known to be below the horizon
     self.navSanityMaxResid    = 25.0                    # meters per SV, normalized nav residuals