GLO acquisition support

Author: Pasi Miettinen

peregrine/acquisition.py

@@ -17,9 +17,11 @@
 import numpy as np
 import pyfftw
 import cPickle
-import defaults
 
 from include.generateCAcode import caCodes
+from include.generateGLOcode import GLOCode
+from peregrine.gps_constants import L1CA
+from peregrine.glo_constants import GLO_L1, glo_l1_step
 
 import logging
 logger = logging.getLogger(__name__)
@@ -191,13 +193,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
@@ -342,24 +346,39 @@ def find_peak(self, freqs, results, interpolation='gaussian'):
     freq_index, cp_samples = np.unravel_index(results.argmax(),
                                               results.shape)
 
-    if freq_index > 1 and freq_index < len(freqs) - 1:
-      delta = self.interpolate(
-          results[freq_index - 1][cp_samples],
-          results[freq_index][cp_samples],
-          results[freq_index + 1][cp_samples],
-          interpolation
-      )
-      if delta > 0:
-        freq = freqs[freq_index] + \
-            (freqs[freq_index + 1] - freqs[freq_index]) * delta
-      else:
-        freq = freqs[freq_index] - \
-            (freqs[freq_index - 1] - freqs[freq_index]) * delta
-    else:
-      freq = freqs[freq_index]
-
     code_phase = float(cp_samples) / self.samples_per_chip
 
+    #print results[freq_index]
+    #pp(results[freq_index - 6][cp_samples])
+    #pp(results[freq_index - 5][cp_samples])
+    #pp(results[freq_index - 4][cp_samples])
+    #pp(results[freq_index - 3][cp_samples])
+    #pp(results[freq_index - 2][cp_samples])
+    #pp(results[freq_index - 1][cp_samples])
+    #pp(results[freq_index][cp_samples])
+    #pp(results[freq_index + 1][cp_samples])
+    #pp(results[freq_index + 2][cp_samples])
+    #pp(results[freq_index + 3][cp_samples])
+    #pp(results[freq_index + 4][cp_samples])
+    #pp(results[freq_index + 5][cp_samples])
+    #pp(results[freq_index + 6][cp_samples])
+
+    #if freq_index > 3 and freq_index < len(freqs) - 3:
+      #delta = self.interpolate(
+          #results[freq_index - 3][cp_samples],
+          #results[freq_index][cp_samples],
+          #results[freq_index + 3][cp_samples],
+          #interpolation
+      #)
+      #if delta > 0:
+        #freq = freqs[freq_index] + \
+            #(freqs[freq_index + 3] - freqs[freq_index]) * delta
+      #else:
+        #freq = freqs[freq_index] - \
+            #(freqs[freq_index - 3] - freqs[freq_index]) * delta
+    #else:
+    freq = freqs[freq_index]
+
     # Calculate SNR for the peak.
     results_mean = np.mean(results)
     if results_mean != 0:
@@ -370,7 +389,8 @@ def find_peak(self, freqs, results, interpolation='gaussian'):
     return (code_phase, freq, snr)
 
   def acquisition(self,
-                  prns=range(32),
+                  prns=xrange(32),
+                  channels=[x - 7 for x in xrange(14)],
                   doppler_priors=None,
                   doppler_search=7000,
                   doppler_step=None,
@@ -379,10 +399,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
@@ -394,8 +414,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.
@@ -413,10 +438,11 @@ 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")
+    logger.info("Acquisition starting for " + self.signal)
     from peregrine.parallel_processing import parmap
 
     # If the Doppler step is not specified, compute it from the coarse
@@ -428,9 +454,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
@@ -446,33 +469,55 @@ def acquisition(self,
       show_progress = False
       logger.warning("show_progress = True but progressbar module not found.")
 
+    if self.signal == 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)),
+                                     maxval=int(input_len *
+                                       (2 * doppler_search / doppler_step + 1)),
                                      fd=progress_fd)
       pbar.start()
     else:
       pbar = None
 
     def do_acq(n):
-      prn = prns[n]
+      if self.signal == L1CA:
+        prn = prns[n]
+        code = caCodes[prn]
+        int_f = self.IF
+        attr = {'prn': prn + 1}
+      else:
+        ch = channels[n]
+        code = GLOCode
+        int_f = self.IF + ch * glo_l1_step
+        attr = {'ch': ch}
       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,
@@ -485,13 +530,22 @@ 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,
-                                     self.signal)
+      if self.signal == L1CA:
+        acq_result = AcquisitionResult(prn,
+                                       carr_freq,
+                                       carr_freq - int_f,
+                                       code_phase,
+                                       snr,
+                                       status,
+                                       L1CA)
+      else:
+        acq_result = GloAcquisitionResult(ch,
+                                          carr_freq,
+                                          carr_freq - int_f,
+                                          code_phase,
+                                          snr,
+                                          status,
+                                          GLO_L1)
 
       # If the acquisition was successful, log it
       if (snr > threshold):
@@ -501,9 +555,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
 
@@ -512,9 +566,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 self.signal == L1CA:
+      logger.info("Acquired %d satellites, PRNs: %s.", len(acq), acq)
+    else:
+      logger.info("Acquired %d channels: %s.", len(acq), acq)
 
     return acq_results
 
@@ -531,7 +587,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.
 
@@ -560,7 +616,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):
@@ -574,7 +630,7 @@ def __init__(self, prn, carr_freq, doppler, code_phase, snr, status, signal,
     self.sample_index = sample_index
 
   def __str__(self):
-    return "PRN %2d (%s) SNR %6.2f @ CP %6.1f, %+8.2f Hz %s" % \
+    return "PRN %2d (%s) SNR %6.2f @ CP %6.3f, %+8.2f Hz %s" % \
         (self.prn + 1, self.signal, self.snr, self.code_phase,
          self.doppler, self.status)
 
@@ -615,6 +671,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.3f, %+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.
@@ -676,4 +746,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

@@ -17,7 +17,7 @@
 ephemMaxAge          = 4 * 3600.0 # Reject an ephemeris entry if older than this
 
 # the size of the sample data block processed at a time
-processing_block_size = 20e6 # [samples]
+processing_block_size = 20e6  # [samples]
 
 # original
 sample_channel_GPS_L1 = 0
@@ -75,18 +75,21 @@
 freq_profile_low_rate = {
     'GPS_L1_IF': 14.58e5,
     'GPS_L2_IF': 7.4e5,
+    'GLO_L1_IF': 12e5,
     'sampling_freq': 24.84375e5}
 
 # 'normal_rate' frequencies profile
 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/glo_constants.py

@@ -0,0 +1,13 @@
+# -*- coding: utf-8 -*-
+from gps_constants import c
+# GLO system parameters
+glo_l1 = 1.602e9  # Hz
+glo_l2 = 1.246e9  # Hz
+glo_code_len = 511
+glo_chip_rate = 0.511e6  # Hz
+glo_l1_step = 0.5625e6  # Hz
+
+glo_code_period = glo_code_len / glo_chip_rate
+glo_code_wavelength = glo_code_period * c
+
+GLO_L1 = 'glo_l1'