diff --git a/.gitignore b/.gitignore index 0b17b6a..0409b13 100644 --- a/.gitignore +++ b/.gitignore @@ -50,3 +50,9 @@ pip-delete-this-directory.txt tests/peregrine_ci_test_data.tar.gz tests/test_data/ tests/test_data_old* + +# Eclipse project files +.pydevproject +.project +.cproject +.settings/ diff --git a/.gitmodules b/.gitmodules index 540aac8..16ea5a7 100644 --- a/.gitmodules +++ b/.gitmodules @@ -1,3 +1,3 @@ [submodule "libswiftnav"] path = libswiftnav - url = https://github.com/swift-nav/libswiftnav + url = https://github.com/adel-mamin/libswiftnav diff --git a/deps.sh b/deps.sh index 33b4f15..37a3d5e 100755 --- a/deps.sh +++ b/deps.sh @@ -6,7 +6,7 @@ function install_deps_ubuntu_maybe () { # Sudo'd version of travis installation instructions sudo apt-get update -qq sudo apt-get install python-software-properties - sudo add-apt-repository --yes ppa:kalakris/cmake + #sudo add-apt-repository --yes ppa:kalakris/cmake sudo apt-get update -qq sudo apt-get -y install cmake \ check \ @@ -17,8 +17,9 @@ function install_deps_ubuntu_maybe () { python-numpy \ python-dev \ cython \ - python-cython \ - python-dev + python-dev \ + python-matplotlib + sudo pip install -U cython git submodule update --init # Build libswiftnav cd libswiftnav/ diff --git a/libswiftnav b/libswiftnav index 2c0d0a7..19a3718 160000 --- a/libswiftnav +++ b/libswiftnav @@ -1 +1 @@ -Subproject commit 2c0d0a776ebdb252ccf94da316c1bf5183b939b0 +Subproject commit 19a3718005fc4c8a8e2fe4a60c381bb5829e34c5 diff --git a/peregrine/acquisition.py b/peregrine/acquisition.py index 01da32d..da0e17e 100644 --- a/peregrine/acquisition.py +++ b/peregrine/acquisition.py @@ -13,12 +13,15 @@ """ +import sys 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, GLO_L1, glo_l1_step import logging logger = logging.getLogger(__name__) @@ -76,12 +79,12 @@ class Acquisition: def __init__(self, samples, - sampling_freq=defaults.sampling_freq, - IF=defaults.IF, - samples_per_code=defaults.samples_per_code, + sampling_freq, + IF, + samples_per_code, code_length=defaults.code_length, n_codes_integrate=4, - offsets = None, + offsets=None, wisdom_file=DEFAULT_WISDOM_FILE): self.sampling_freq = sampling_freq @@ -95,11 +98,11 @@ def __init__(self, if n_codes_integrate <= 10: offsets = [0, self.n_integrate] elif n_codes_integrate <= 13: - offsets = [0, 2*(n_codes_integrate - 10)*self.samples_per_code, + offsets = [0, 2 * (n_codes_integrate - 10) * self.samples_per_code, self.n_integrate] elif n_codes_integrate <= 15: offsets = [0, (n_codes_integrate - 10) * self.samples_per_code, - 2*(n_codes_integrate - 10) * self.samples_per_code, + 2 * (n_codes_integrate - 10) * self.samples_per_code, self.n_integrate] else: raise ValueError("Integration interval too long to guess nav-declobber " @@ -129,9 +132,9 @@ def __init__(self, # Allocate aligned arrays for the inverse FFT. self.corr_ft = pyfftw.n_byte_align_empty((self.n_integrate), 16, - dtype=np.complex128) + dtype=np.complex128) self.corr = pyfftw.n_byte_align_empty((self.n_integrate), 16, - dtype=np.complex128) + dtype=np.complex128) # Setup FFTW transforms for inverse FFT. self.corr_ifft = pyfftw.FFTW(self.corr_ft, self.corr, @@ -187,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 @@ -226,13 +231,13 @@ def interpolate(self, S_0, S_1, S_2, interpolation='gaussian'): """ if interpolation == 'parabolic': # Parabolic interpolation. - return 0.5 * (S_2 - S_0) / (2*S_1 - S_0 - S_2) + return 0.5 * (S_2 - S_0) / (2 * S_1 - S_0 - S_2) elif interpolation == 'gaussian': # Gaussian interpolation. ln_S_0 = np.log(S_0) ln_S_1 = np.log(S_1) ln_S_2 = np.log(S_2) - return 0.5 * (ln_S_2 - ln_S_0) / (2*ln_S_1 - ln_S_0 - ln_S_2) + return 0.5 * (ln_S_2 - ln_S_0) / (2 * ln_S_1 - ln_S_0 - ln_S_2) elif interpolation == 'none': return 0 else: @@ -270,8 +275,9 @@ def acquire(self, code, freqs, progress_callback=None): # Upsample the code to our sampling frequency. code_indices = np.arange(1.0, self.n_integrate + 1.0) / \ - self.samples_per_chip - code_indices = np.remainder(np.asarray(code_indices, np.int), self.code_length) + self.samples_per_chip + code_indices = np.remainder( + np.asarray(code_indices, np.int), self.code_length) self.code[:] = code[code_indices] # Find the conjugate Fourier transform of the code which will be used to @@ -288,7 +294,7 @@ def acquire(self, code, freqs, progress_callback=None): # Shift the signal in the frequency domain to remove the carrier # i.e. mix down to baseband. shift = int((float(freq) * len(self.short_samples_ft[0]) / - self.sampling_freq) + 0.5) + self.sampling_freq) + 0.5) # Search over the possible nav bit offset intervals for offset_i in range(len(self.offsets)): @@ -307,7 +313,6 @@ def acquire(self, code, freqs, progress_callback=None): max_indices = np.unravel_index(results.argmax(), results.shape) return results[max_indices[0]] - def find_peak(self, freqs, results, interpolation='gaussian'): """ Find the peak within an set of acquisition results. @@ -338,17 +343,19 @@ 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: + 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 + 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 + 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 + freq = freqs[freq_index] - \ + (freqs[freq_index - 1] - freqs[freq_index]) * delta else: freq = freqs[freq_index] @@ -360,19 +367,21 @@ def find_peak(self, freqs, results, interpolation='gaussian'): return (code_phase, freq, snr) def acquisition(self, - prns=range(32), - doppler_priors = None, - doppler_search = 7000, - doppler_step = None, + bandcode=L1CA, + prns=xrange(32), + channels=[x - 7 for x in xrange(14)], + doppler_priors=None, + doppler_search=7000, + doppler_step=None, threshold=DEFAULT_THRESHOLD, - show_progress=True, + progress_bar_output='none', 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 @@ -384,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. @@ -403,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 @@ -418,46 +437,74 @@ 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 + elif progress_bar_output == 'stderr': + show_progress = True + progress_fd = sys.stderr + else: + show_progress = False + progress_fd = -1 # If progressbar is not available, disable show_progress. if show_progress and not _progressbar_available: 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 = GLOCode + 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, - interpolation = 'gaussian') + interpolation='gaussian') # If the result is above the threshold, then we have acquired the # satellite. @@ -466,12 +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) + 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): @@ -480,9 +538,10 @@ def progress_callback(freq_num, num_freqs): return acq_result if multi: - acq_results = parmap(do_acq, range(len(prns)), show_progress=show_progress) + acq_results = parmap( + 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 @@ -491,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 @@ -506,10 +567,11 @@ def load_wisdom(self, wisdom_file=DEFAULT_WISDOM_FILE): def save_wisdom(self, wisdom_file=DEFAULT_WISDOM_FILE): """Save FFTW wisdom to file.""" with open(wisdom_file, 'wb') as f: - cPickle.dump(pyfftw.export_wisdom(), f, protocol=cPickle.HIGHEST_PROTOCOL) + cPickle.dump( + pyfftw.export_wisdom(), f, protocol=cPickle.HIGHEST_PROTOCOL) -class AcquisitionResult: +class AcquisitionResult(object): """ Stores the acquisition parameters of a single satellite. @@ -531,22 +593,30 @@ class AcquisitionResult: * `'A'` : The satellite has been successfully acquired. * `'-'` : The acquisition was not successful, the SNR was below the acquisition threshold. - + signal : {'l1ca', 'l2c'} + The type of the signal: L1C/A or L2C + sample_channel : IQ channel index + sample_index : Index of sample when acquisition succeeded """ - __slots__ = ('prn', 'carr_freq', 'doppler', 'code_phase', 'snr', 'status') + __slots__ = ('prn', 'carr_freq', 'doppler', + 'code_phase', 'snr', 'status', 'signal', 'sample_index') - def __init__(self, prn, carr_freq, doppler, code_phase, snr, status): + def __init__(self, prn, carr_freq, doppler, code_phase, snr, status, signal, + sample_index=0): self.prn = prn self.snr = snr self.carr_freq = carr_freq self.doppler = doppler self.code_phase = code_phase self.status = status + self.signal = signal + self.sample_index = sample_index def __str__(self): - return "PRN %2d SNR %6.2f @ CP %6.1f, %+8.2f Hz %s" % \ - (self.prn + 1, self.snr, self.code_phase, self.doppler, self.status) + return "PRN %2d (%s) SNR %6.2f @ CP %6.1f, %+8.2f Hz %s" % \ + (self.prn + 1, self.signal, self.snr, self.code_phase, + self.doppler, self.status) def __repr__(self): return "" % self.__str__() @@ -570,7 +640,7 @@ def _equal(self, other): ------ out : bool True if the passed :class:`AcquisitionResult` object is identical. - + """ if set(self.__dict__.keys()) != set(other.__dict__.keys()): return False @@ -584,6 +654,21 @@ 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. @@ -599,6 +684,7 @@ def save_acq_results(filename, acq_results): with open(filename, 'wb') as f: cPickle.dump(acq_results, f, protocol=cPickle.HIGHEST_PROTOCOL) + def load_acq_results(filename): """ Load a set of acquisition results from a file. @@ -617,7 +703,8 @@ def load_acq_results(filename): with open(filename, 'rb') as f: return cPickle.load(f) -def print_scores(acq_results, pred, pred_dopp = None): + +def print_scores(acq_results, pred, pred_dopp=None): if pred_dopp is None: pred_dopp = np.zeros_like(pred) @@ -628,19 +715,20 @@ def print_scores(acq_results, pred, pred_dopp = None): sum_abs_dopp_err = 0 for i, prn in enumerate(pred): - print "%2d\t%+6.0f" % (prn + 1, pred_dopp[i]), - if acq_results[i].status == 'A': - n_match += 1 - dopp_err = acq_results[i].doppler - pred_dopp[i] - sum_dopp_err += dopp_err - sum_abs_dopp_err += abs(dopp_err) - if abs(dopp_err) > abs(worst_dopp_err): - worst_dopp_err = dopp_err - print "\t%+6.0f\t%+5.0f\t%5.1f" % ( - acq_results[i].doppler, dopp_err, acq_results[i].snr) - else: - print + print "%2d\t%+6.0f" % (prn + 1, pred_dopp[i]), + if acq_results[i].status == 'A': + n_match += 1 + dopp_err = acq_results[i].doppler - pred_dopp[i] + sum_dopp_err += dopp_err + sum_abs_dopp_err += abs(dopp_err) + if abs(dopp_err) > abs(worst_dopp_err): + worst_dopp_err = dopp_err + print "\t%+6.0f\t%+5.0f\t%5.1f" % ( + acq_results[i].doppler, dopp_err, acq_results[i].snr) + else: + print 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) diff --git a/peregrine/analysis/print_res.py b/peregrine/analysis/print_res.py new file mode 100755 index 0000000..48c74b9 --- /dev/null +++ b/peregrine/analysis/print_res.py @@ -0,0 +1,53 @@ +#!/usr/bin/env python + +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Adel Mamin +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + +from scipy import signal +import numpy as np +import matplotlib.pyplot as plt +from StringIO import StringIO +import argparse + + +def main(): + parser = argparse.ArgumentParser() + + parser.add_argument("-f", "--file", default="tracking_res.csv", + help="the input CSV file to process") + + parser.add_argument("-p", "--par-to-print", default="CN0", + help="parameter to print") + + parser.add_argument("-s", "--time-step", default="0.1", + help="time step [s]") + + args = parser.parse_args() + + fig = plt.figure() + #plt.title('Carrier tracking loop filter frequency response') + ax1 = fig.add_subplot(111) + + plt.ylabel(args.par_to_print, color='b') + plt.xlabel('time [%s s]' % float(args.time_step)) + + data = np.genfromtxt(args.file, dtype=float, delimiter=',', names=True) + + plt.plot(np.array(range(len(data[args.par_to_print]))), + np.array(data[args.par_to_print]), 'r.') + + plt.legend(loc='upper right') + + plt.grid() + plt.axis('tight') + plt.show() + +if __name__ == '__main__': + main() diff --git a/peregrine/analysis/samples.py b/peregrine/analysis/samples.py index efeb553..0ca1821 100755 --- a/peregrine/analysis/samples.py +++ b/peregrine/analysis/samples.py @@ -59,7 +59,7 @@ def hist(samples, ax=None, value_range=None, bin_width=1.0, max_len=ANALYSIS_MAX """ if max_len is not None and len(samples) > max_len: - logger.debug( "Truncating to %d samples." % max_len) + logger.debug("Truncating to %d samples." % max_len) samples = samples[:max_len] if ax is None: @@ -72,8 +72,8 @@ def hist(samples, ax=None, value_range=None, bin_width=1.0, max_len=ANALYSIS_MAX max_val = np.max(samples) min_val = np.min(samples) - n_bins = 1 + np.round(float(max_val) - float(min_val) / bin_width) - bins = np.linspace(min_val-bin_width/2.0, max_val+bin_width/2.0, n_bins+1) + n_bins = 1 + np.round(float(max_val) - float(min_val) / bin_width) + bins = np.linspace(min_val - bin_width / 2.0, max_val + bin_width / 2.0, n_bins + 1) ticks = np.linspace(min_val, max_val, n_bins) @@ -83,11 +83,11 @@ def hist(samples, ax=None, value_range=None, bin_width=1.0, max_len=ANALYSIS_MAX ax.set_xlabel('Sample value') if len(ticks) < 22: ax.set_xticks(ticks) - ax.set_xbound(min_val-bin_width, max_val+bin_width) + ax.set_xbound(min_val - bin_width, max_val + bin_width) ax.set_ylabel('Count') y_min, y_max = ax.get_ybound() - ax.set_ybound(0, y_max*1.1) - ax.ticklabel_format(style='sci', scilimits=(0,0), axis='y') + ax.set_ybound(0, y_max * 1.1) + ax.ticklabel_format(style='sci', scilimits=(0, 0), axis='y') return ax @@ -118,7 +118,7 @@ def psd(samples, sampling_freq=None, ax=None, max_len=ANALYSIS_MAX_LEN): """ if max_len is not None and len(samples) > max_len: - logger.debug( "Truncating to %d samples." % max_len) + logger.debug("Truncating to %d samples." % max_len) samples = samples[:max_len] if ax is None: @@ -132,7 +132,7 @@ def psd(samples, sampling_freq=None, ax=None, max_len=ANALYSIS_MAX_LEN): ax.set_ylabel('Power Spectral Density ($f_s \cdot \mathrm{Hz}^{-1}$)') sampling_freq = 1.0 else: - ax.ticklabel_format(style='sci', scilimits=(0,0), axis='x') + ax.ticklabel_format(style='sci', scilimits=(0, 0), axis='x') ax.set_xlabel('Frequency (Hz)') ax.set_ylabel('Power Spectral Density ($\mathrm{Hz}^{-1}$)') @@ -143,7 +143,7 @@ def psd(samples, sampling_freq=None, ax=None, max_len=ANALYSIS_MAX_LEN): noverlap=1024) ax.semilogy(freqs, Pxx, color='black') - ax.set_xbound(0, sampling_freq/2.0) + ax.set_xbound(0, sampling_freq / 2.0) return ax @@ -172,8 +172,8 @@ def summary(samples, sampling_freq=None, max_len=ANALYSIS_MAX_LEN): ax1 = fig.add_subplot(121) ax2 = fig.add_subplot(122) - hist(samples, ax=ax1, max_len=max_len) - psd(samples, sampling_freq, ax=ax2, max_len=max_len) + hist(samples[0], ax=ax1, max_len=max_len) + psd(samples[0], sampling_freq, ax=ax2, max_len=max_len) fig.set_size_inches(10, 4, forward=True) fig.tight_layout() @@ -193,7 +193,7 @@ def main(): + "'int8', '1bit', '1bitrev' or 'piksi' (default)") args = parser.parse_args() - samples = peregrine.samples.load_samples(args.file, args.num_samples, file_format = args.format) + samples = peregrine.samples.load_samples(args.file, args.num_samples, file_format=args.format) summary(samples) plt.show() diff --git a/peregrine/analysis/tracking_loop.py b/peregrine/analysis/tracking_loop.py new file mode 100755 index 0000000..51fd650 --- /dev/null +++ b/peregrine/analysis/tracking_loop.py @@ -0,0 +1,270 @@ +#!/usr/bin/env python + +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Adel Mamin +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + +import argparse +import sys +from peregrine.samples import load_samples +from peregrine.acquisition import AcquisitionResult +from peregrine import defaults +from peregrine.log import default_logging_config +from peregrine.tracking import Tracker +from peregrine.gps_constants import L1CA, L2C +from peregrine.initSettings import initSettings + + +def main(): + default_logging_config() + + parser = argparse.ArgumentParser() + + if sys.stdout.isatty(): + progress_bar_default = 'stdout' + elif sys.stderr.isatty(): + progress_bar_default = 'stderr' + else: + progress_bar_default = 'none' + + parser.add_argument("--progress-bar", + metavar='FLAG', + choices=['stdout', 'stderr', 'none'], + default=progress_bar_default, + help="Show progress bar. Default is '%s'" % + progress_bar_default) + + inputCtrl = parser.add_argument_group('Data Source', + 'Data source configuration') + inputCtrl.add_argument("file", + help="The sample data file to process") + + inputCtrl.add_argument("-f", "--file-format", + choices=['piksi', 'int8', '1bit', '1bitrev', + '1bit_x2', '2bits', '2bits_x2', '2bits_x4'], + metavar='FORMAT', + help="The format of the sample data file " + "('piksi', 'int8', '1bit', '1bitrev', " + "'1bit_x2', '2bits', '2bits_x2', '2bits_x4')") + + inputCtrl.add_argument("-t", "--ms-to-track", + metavar='MS', + help="the number of milliseconds to track." + "(-1: use all available data", + default="-1") + + inputCtrl.add_argument("--skip-samples", + default=0, + metavar='N_SAMPLES', + help="How many samples to skip") + + inputCtrl.add_argument("-s", + "--sampling-freq", + metavar='FREQ', + help="Sampling frequency [Hz]. ") + + inputCtrl.add_argument("--profile", + choices=['peregrine', 'custom_rate', 'low_rate', + 'normal_rate', 'piksi_v3', 'high_rate'], + metavar='PROFILE', + help="L1C/A & L2C IF + sampling frequency profile" + "('peregrine'/'custom_rate', 'low_rate', " + "'normal_rate', 'piksi_v3', 'high_rate')", + default='peregrine') + + signalParam = parser.add_argument_group('Signal tracking', + 'Parameters for satellite vehicle' + ' signal') + + signalParam.add_argument("-P", "--prn", + help="PRN to track. ") + + signalParam.add_argument("-p", "--code-phase", + metavar='CHIPS', + help="Code phase [chips]. ") + + signalParam.add_argument("-d", "--carr-doppler", + metavar='DOPPLER', + help="carrier Doppler frequency [Hz]. ") + + signalParam.add_argument("-S", "--signal", + choices=[L1CA, L2C], + metavar='BAND', + help="Signal type (l1ca / l2c)") + signalParam.add_argument("--l2c-handover", + action='store_true', + help="Perform L2C handover", + default=False) + signalParam.add_argument('--l1ca-profile', + metavar='PROFILE', + help='L1 C/A stage profile. Controls coherent' + ' integration time and tuning parameters: %s.' % + str(defaults.l1ca_stage_profiles.keys()), + choices=defaults.l1ca_stage_profiles.keys()) + + fpgaSim = parser.add_argument_group('FPGA simulation', + 'FPGA delay control simulation') + fpgaExcl = fpgaSim.add_mutually_exclusive_group(required=False) + fpgaExcl.add_argument("--pipelining", + type=float, + nargs='?', + metavar='PIPELINING_K', + help="Use FPGA pipelining simulation. Supply optional " + " coefficient (%f)" % defaults.pipelining_k, + const=defaults.pipelining_k, + default=None) + + fpgaExcl.add_argument("--short-long-cycles", + type=float, + nargs='?', + metavar='PIPELINING_K', + help="Use FPGA short-long cycle simulation. Supply" + " optional pipelining coefficient (0.)", + const=0., + default=None) + + outputCtrl = parser.add_argument_group('Output parameters', + 'Parameters that control output' + ' data stream.') + outputCtrl.add_argument("-o", "--output-file", + default="track.csv", + help="Track results file name. Default: %s" % + "track.csv") + + args = parser.parse_args() + + if args.file is None: + parser.print_help() + return + + skip_samples = int(args.skip_samples) + + if args.profile == 'peregrine' or args.profile == 'custom_rate': + freq_profile = defaults.freq_profile_peregrine + elif args.profile == 'low_rate': + freq_profile = defaults.freq_profile_low_rate + elif args.profile == 'normal_rate': + freq_profile = defaults.freq_profile_normal_rate + elif args.profile == 'high_rate': + freq_profile = defaults.freq_profile_high_rate + else: + raise NotImplementedError() + + isL1CA = (args.signal == L1CA) + isL2C = (args.signal == L2C) + + if isL1CA: + signal = L1CA + IF = freq_profile['GPS_L1_IF'] + elif isL2C: + signal = L2C + IF = freq_profile['GPS_L2_IF'] + else: + raise NotImplementedError() + + if args.l2c_handover and not isL2C: + l2c_handover = True + else: + l2c_handover = False + + if args.sampling_freq is not None: + sampling_freq = float(args.sampling_freq) # [Hz] + else: + sampling_freq = freq_profile['sampling_freq'] # [Hz] + + # Initialize constants, settings + settings = initSettings(freq_profile) + + settings.fileName = args.file + + carr_doppler = float(args.carr_doppler) + code_phase = float(args.code_phase) + prn = int(args.prn) - 1 + + ms_to_track = int(args.ms_to_track) + + if args.pipelining is not None: + tracker_options = {'mode': 'pipelining', + 'k': args.pipelining} + elif args.short_long_cycles is not None: + tracker_options = {'mode': 'short-long-cycles', + 'k': args.short_long_cycles} + else: + tracker_options = None + + acq_result = AcquisitionResult(prn=prn, + snr=25, # dB + carr_freq=IF + carr_doppler, + doppler=carr_doppler, + code_phase=code_phase, + status='A', + signal=signal, + sample_index=skip_samples) + + if args.l1ca_profile: + profile = defaults.l1ca_stage_profiles[args.l1ca_profile] + stage2_coherent_ms = profile[1]['coherent_ms'] + stage2_params = profile[1]['loop_filter_params'] + else: + stage2_coherent_ms = None + stage2_params = None + + samples = {L1CA: {'IF': freq_profile['GPS_L1_IF']}, + L2C: {'IF': freq_profile['GPS_L2_IF']}, + 'samples_total': -1, + 'sample_index': skip_samples} + + load_samples(samples=samples, + filename=args.file, + file_format=args.file_format) + + if ms_to_track < 0: + # use all available data + ms_to_track = int(1e3 * samples['samples_total'] / sampling_freq) + + print "==================== Tracking parameters =============================" + print "File: %s" % args.file + print "File format: %s" % args.file_format + print "PRN to track [1-32]: %s" % args.prn + print "Time to process [ms]: %s" % ms_to_track + print "L1 IF [Hz]: %f" % freq_profile['GPS_L1_IF'] + print "L2 IF [Hz]: %f" % freq_profile['GPS_L2_IF'] + print "Sampling frequency [Hz]: %f" % sampling_freq + print "Initial carrier Doppler frequency [Hz]: %s" % carr_doppler + print "Initial code phase [chips]: %s" % code_phase + print "Signal: %s" % args.signal + print "L1 stage profile: %s" % args.l1ca_profile + print "Tracker options: %s" % str(tracker_options) + print "======================================================================" + + tracker = Tracker(samples=samples, + channels=[acq_result], + ms_to_track=ms_to_track, + sampling_freq=sampling_freq, # [Hz] + l2c_handover=l2c_handover, + stage2_coherent_ms=stage2_coherent_ms, + stage2_loop_filter_params=stage2_params, + tracker_options=tracker_options, + output_file=args.output_file, + progress_bar_output=args.progress_bar) + tracker.start() + condition = True + while condition: + sample_index = tracker.run_channels(samples) + if sample_index == samples['sample_index']: + condition = False + else: + samples['sample_index'] = sample_index + load_samples(samples=samples, + filename=args.file, + file_format=args.file_format) + tracker.stop() + +if __name__ == '__main__': + main() diff --git a/peregrine/defaults.py b/peregrine/defaults.py index 8fccec7..8e9e35c 100644 --- a/peregrine/defaults.py +++ b/peregrine/defaults.py @@ -1,4 +1,5 @@ # Copyright (C) 2014 Swift Navigation Inc. +# Contact: Adel Mamin # # This source is subject to the license found in the file 'LICENSE' which must # be be distributed together with this source. All other rights reserved. @@ -7,14 +8,255 @@ # EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED # WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. -ms_to_track = 37*1e3 +ms_to_track = 37 * 1e3 skip_samples = 1000 file_format = 'piksi' -IF = 4.092e6 # Hz -sampling_freq = 16.368e6 # Hz -chipping_rate = 1.023e6 # Hz -code_length = 1023 # chips +processing_block_size = 20e6 # [samples] + +chipping_rate = 1.023e6 # Hz +code_length = 1023 # chips code_period = code_length / chipping_rate -samples_per_code = code_period * sampling_freq +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 +sample_channel_GLO_L1 = 2 +sample_channel_GLO_L2 = 3 + +file_encoding_1bit_x2 = [ + sample_channel_GPS_L1, # GPS L1 + sample_channel_GPS_L2] # GPS L2 + +file_encoding_2bits_x2 = file_encoding_1bit_x2 + +# encoding is taken from here: +# https://swiftnav.hackpad.com/MicroZed-Sample-Grabber-IFgt5DbAunD +# 2 bits per frontend channel in every byte: +# RF4 RF3 RF2 RF1 +# 00 00 00 00 +# +# RF 1: +# GPS L1 @ 14.58MHz (1575.42MHz) +# Galileo E1 @ 14.58MHz (1575.42MHz) +# Beidou B1 @ 28.902 MHz (1561.098 MHz) +# +# RF 2: +# GLONASS L1 @ 12MHz (1602MHz) +# +# RF 3: +# GLONASS L2 @ 11MHz (1246MHz) +# Beidou B3 @ 33.52MHz (1268.52MHz) +# Galileo E6 @ 43.75 MHz(1278.75MHz) +# +# RF 4: +# GPS L2 @ 7.4MHz (1227.6MHz) +# Galileo E5b-I/Q @ 27.86MHz (1207.14MHz) +# Beidou B2 @ 27.86MHz (1207.14MHz) +file_encoding_2bits_x4 = [ + sample_channel_GPS_L2, # RF4 + sample_channel_GLO_L2, # RF3 + sample_channel_GLO_L1, # RF2 + sample_channel_GPS_L1] # RF1 + +file_encoding_profile = { + '1bit_x2': file_encoding_1bit_x2, + '2bits_x2': file_encoding_2bits_x2, + '2bits_x4': file_encoding_2bits_x4} + +# 'peregrine' frequencies profile +freq_profile_peregrine = { + 'GPS_L1_IF': 4.092e6, + 'GPS_L2_IF': 4.092e6, + 'sampling_freq': 16.368e6} + +# 'low_rate' frequencies profile +freq_profile_low_rate = { + 'GPS_L1_IF': 14.58e5, + 'GPS_L2_IF': 7.4e5, + '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} + +# '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 +L2C_CHANNEL_BANDWIDTH_HZ = 1000 + +l1ca_stage1_loop_filter_params = { + "loop_freq": 1e3, # loop frequency [Hz] + "code_bw": 1, # Code loop NBW + "code_zeta": 0.7, # Code loop zeta + "code_k": 1, # Code loop k + "carr_to_code": 1540, # Carrier-to-code freq ratio (carrier aiding) + "carr_bw": 10, # Carrier loop NBW + "carr_zeta": 0.7, # Carrier loop zeta + "carr_k": 1, # Carrier loop k + "carr_freq_b1": 5} # Carrier loop aiding_igain + +l2c_loop_filter_params = { + "loop_freq": 50, # loop frequency [Hz] + "code_bw": 1.4, # Code loop NBW + "code_zeta": 0.707, # Code loop zeta + "code_k": 1, # Code loop k + "carr_to_code": 1200, # Carrier-to-code freq ratio (carrier aiding) + "carr_bw": 13, # Carrier loop NBW + "carr_zeta": 0.707, # Carrier loop zeta + "carr_k": 1, # Carrier loop k + "carr_freq_b1": 5} # Carrier loop aiding_igain + + +# Tracking stages. See track.c for more details. +# 1;20 ms stages +l1ca_stage_params_slow = \ + ({'coherent_ms': 1, + 'loop_filter_params': {'code_params': (1., 0.7, 1.), # NBW, zeta, k + 'carr_params': (10., 0.7, 1.), # NBW, zeta, k + 'loop_freq': 1000., # 1000/coherent_ms + 'carr_freq_igain': 5., # fll_aid + 'carr_to_code': 1540. # carr_to_code + } + }, + {'coherent_ms': 20, + 'loop_filter_params': {'code_params': (1., 0.7, 1.), # NBW, zeta, k + 'carr_params': (12., 0.7, 1.), # NBW, zeta, k + 'loop_freq': 1000. / 20, # 1000/coherent_ms + 'carr_freq_igain': 0., # fll_aid + 'carr_to_code': 1540. # carr_to_code + } + } + ) + +l1ca_stage_params_slow2 = \ + ({'coherent_ms': 1, + 'loop_filter_params': {'code_params': (1., 0.7, 1.), # NBW, zeta, k + 'carr_params': (10., 0.7, 1.), # NBW, zeta, k + 'loop_freq': 1000., # 1000/coherent_ms + 'carr_freq_igain': 5., # fll_aid + 'carr_to_code': 1540. # carr_to_code + } + }, + {'coherent_ms': 10, + 'loop_filter_params': {'code_params': (1., 0.7, 1.), # NBW, zeta, k + 'carr_params': (12., 0.7, 1.), # NBW, zeta, k + 'loop_freq': 1000. / 10, # 1000/coherent_ms + 'carr_freq_igain': 0., # fll_aid + 'carr_to_code': 1540. # carr_to_code + } + } + ) + +# 1;5 ms stages +l1ca_stage_params_med = \ + ({'coherent_ms': 1, + 'loop_filter_params': {'code_params': (1., 0.7, 1.), # NBW, zeta, k + 'carr_params': (10., 0.7, 1.), # NBW, zeta, k + 'loop_freq': 1000., # 1000/coherent_ms + 'carr_freq_igain': 5., # fll_aid + 'carr_to_code': 1540. # carr_to_code + } + }, + + {'coherent_ms': 5, + 'loop_filter_params': {'code_params': (1., 0.7, 1.), # NBW, zeta, k + 'carr_params': (50., 0.7, 1.), # NBW, zeta, k + 'loop_freq': 1000. / 5, # 1000/coherent_ms + 'carr_freq_igain': 0., # fll_aid + 'carr_to_code': 1540. # carr_to_code + } + } + ) + +# 1;4 ms stages +l1ca_stage_params_fast = \ + ({'coherent_ms': 1, + 'loop_filter_params': {'code_params': (1., 0.7, 1.), # NBW, zeta, k + 'carr_params': (10., 0.7, 1.), # NBW, zeta, k + 'loop_freq': 1000., # 1000/coherent_ms + 'carr_freq_igain': 5., # fll_aid + 'carr_to_code': 1540. # carr_to_code + } + }, + {'coherent_ms': 4, + 'loop_filter_params': {'code_params': (1., 0.7, 1.), # NBW, zeta, k + 'carr_params': (62., 0.7, 1.), # NBW, zeta, k + 'loop_freq': 1000. / 4, # 1000/coherent_ms + 'carr_freq_igain': 0., # fll_aid + 'carr_to_code': 1540. # carr_to_code + } + } + ) + +# 1;2 ms stages +l1ca_stage_params_extrafast = \ + ({'coherent_ms': 1, + 'loop_filter_params': {'code_params': (1., 0.7, 1.), # NBW, zeta, k + 'carr_params': (10., 0.7, 1.), # NBW, zeta, k + 'loop_freq': 1000., # 1000/coherent_ms + 'carr_freq_igain': 5., # fll_aid + 'carr_to_code': 1540. # carr_to_code + } + }, + {'coherent_ms': 2, + 'loop_filter_params': {'code_params': (1., 0.7, 1.), # NBW, zeta, k + 'carr_params': (100., 0.7, 1.), # NBW, zeta, k + 'loop_freq': 1000. / 2, # 1000/coherent_ms + 'carr_freq_igain': 0., # fll_aid + 'carr_to_code': 1540. # carr_to_code + } + } + ) + +# L1 C/A stage profiles +l1ca_stage_profiles = {'slow': l1ca_stage_params_slow, + 'slow2': l1ca_stage_params_slow2, + 'med': l1ca_stage_params_med, + 'fast': l1ca_stage_params_fast, + 'extrafast': l1ca_stage_params_extrafast} + +# pessimistic set +l1ca_lock_detect_params_pess = {"k1": 0.10, "k2": 1.4, "lp": 200, "lo": 50} + +# normal set +l1ca_lock_detect_params_normal = {"k1": 0.05, "k2": 1.4, "lp": 150, "lo": 50} + +# optimal set +l1ca_lock_detect_params_opt = {"k1": 0.02, "k2": 1.1, "lp": 150, "lo": 50} + +# extra optimal set +l1ca_lock_detect_params_extraopt = {"k1": 0.02, "k2": 0.8, "lp": 150, "lo": 50} + +# disable lock detect +l1ca_lock_detect_params_disable = {"k1": 0.02, "k2": 1e-6, "lp": 1, "lo": 1} + + +# L2C 20ms lock detect profile +# References: +# - Understanding GPS: Principles and Applications. +# Elliott D. Kaplan. Artech House, 2006. 2nd edition +# p.235 +l2c_lock_detect_params_20ms = { + 'k1': 0.0247, # LPF with -3dB at ~0.4 Hz + 'k2': 1.5, # use ~26 degrees I/Q phase angle as a threshold + 'lp': 50, # 1000ms worth of I/Q samples to reach pessimistic lock + 'lo': 240} # 4800ms worth of I/Q samples to lower optimistic lock + +alias_detect_interval_ms = 500 + +# Default pipelining prediction coefficient +pipelining_k = .9549 diff --git a/peregrine/gps_constants.py b/peregrine/gps_constants.py index 789a89e..620cc27 100644 --- a/peregrine/gps_constants.py +++ b/peregrine/gps_constants.py @@ -2,19 +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 +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' diff --git a/peregrine/include/generateL2CMcode.py b/peregrine/include/generateL2CMcode.py new file mode 100644 index 0000000..ad83dbd --- /dev/null +++ b/peregrine/include/generateL2CMcode.py @@ -0,0 +1,99 @@ +# Copyright (C) 2016 Swift Navigation Inc. +# +# Contact: Dmitry Tatarinov +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + +import numpy as np + +def generateL2CMcode(PRN): + """ + The function generates PRN sequence for a particular SV. + In the sequence '0' is represented as '1', 1 as -1 + + INPUT: SV number from 0 (SV1) to 31 (SV32) + OUTPUT: + - PRN seequence array of 10230, + - end state of shift register for testing purpuses + """ + #--- Sanity sheck for PRN number ------------------------------------------ + if PRN < 0 or PRN > 31: + raise ValueError('PRN number(',PRN,') is not in range [0..31]') + + #--- Initial states for shift register for each PRN[1..32], --------------- + # see IS-GPS-200H, Table 3-IIa + initL2CM = [\ + 0742417664, #PRN 1 + 0756014035, + 0002747144, + 0066265724, + 0601403471, + 0703232733, + 0124510070, + 0617316361, + 0047541621, + 0733031046, + 0713512145, + 0024437606, + 0021264003, + 0230655351, + 0001314400, + 0222021506, + 0540264026, + 0205521705, + 0064022144, + 0120161274, + 0044023533, + 0724744327, + 0045743577, + 0741201660, + 0700274134, + 0010247261, + 0713433445, + 0737324162, + 0311627434, + 0710452007, + 0722462133, + 0050172213 #PRN 32 + ] + + #--- Init L2CM PRN and shift reg ------------------------------------------ + L2CM_PRN = np.zeros(10230, np.int8) + + #--- Load Shift register -------------------------------------------------- + shift_cm = initL2CM[PRN] + + #--- Generate L2CM PRN chips ---------------------------------------------- + for i in range(10230): + + out = shift_cm & 1 + if out == 1: + L2CM_PRN[i] = -1 #-1 to represent '1' + else: + L2CM_PRN[i] = 1 # 1 to represent '0' + + shift_reg_out = shift_cm + + shift_cm ^= out << 3 + shift_cm ^= out << 4 + shift_cm ^= out << 5 + shift_cm ^= out << 6 + shift_cm ^= out << 9 + shift_cm ^= out << 11 + shift_cm ^= out << 13 + shift_cm ^= out << 16 + shift_cm ^= out << 19 + shift_cm ^= out << 21 + shift_cm ^= out << 24 + shift_cm = (shift_cm >> 1) | (out << 26) + + return (L2CM_PRN, shift_reg_out) + +L2CMCodes = np.empty((32,10230), dtype=np.int8) +for PRN in range(32): + L2CMCodes[PRN][:] = generateL2CMcode(PRN)[0] + diff --git a/peregrine/initSettings.py b/peregrine/initSettings.py index c17cf16..1333633 100644 --- a/peregrine/initSettings.py +++ b/peregrine/initSettings.py @@ -10,17 +10,20 @@ import defaults class initSettings: - def __init__(self): - + def __init__(self, freq_profile): self.msToProcess = 39000 # Number of ms of samples to perform tracking over (ms) self.skipNumberOfBytes = 0 # Skip bytes in sample file before loading samples for acquisition (bytes) - self.IF = defaults.IF # Intermediate frequency of signal in sample file (Hz) - self.samplingFreq = defaults.sampling_freq # Sampling frequency of sample file (Hz) + 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 (meters) + self.navSanityMaxResid = 25.0 # meters per SV, normalized nav residuals self.abortIfInsane = True # Abort the whole attempt if sanity check fails self.useCache = True self.cacheDir = 'cache' diff --git a/peregrine/iqgen/__init__.py b/peregrine/iqgen/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/peregrine/iqgen/bits/__init__.py b/peregrine/iqgen/bits/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/peregrine/iqgen/bits/amplitude_base.py b/peregrine/iqgen/bits/amplitude_base.py new file mode 100644 index 0000000..67a070d --- /dev/null +++ b/peregrine/iqgen/bits/amplitude_base.py @@ -0,0 +1,57 @@ +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Valeri Atamaniouk +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + +""" +The :mod:`peregrine.iqgen.amplitude_base` module contains classes and functions +related to base implementation of amplitude class. + +""" + + +class AmplitudeBase(object): + ''' + Amplitude control for a signal source. + ''' + + def __init__(self): + ''' + Constructs base object for amplitude control. + ''' + super(AmplitudeBase, self).__init__() + + def applyAmplitude(self, signal, userTimeAll_s): + ''' + Applies amplitude modulation to signal. + + Parameters + ---------- + signal : numpy.ndarray + Signal sample vector. Each element defines signal amplitude in range + [-1; +1]. This vector is modified in place. + userTimeAll_s : numpy.ndarray + Sample time vector. Each element defines sample time in seconds. + + Returns + ------- + numpy.ndarray + Array with output samples + ''' + raise NotImplementedError() + + def computeMeanPower(self): + ''' + Computes mean signal power. + + Returns + ------- + float + Mean signal power for the configured amplitude + ''' + raise NotImplementedError() diff --git a/peregrine/iqgen/bits/amplitude_factory.py b/peregrine/iqgen/bits/amplitude_factory.py new file mode 100644 index 0000000..1c88ea3 --- /dev/null +++ b/peregrine/iqgen/bits/amplitude_factory.py @@ -0,0 +1,68 @@ +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Valeri Atamaniouk +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + +""" +The :mod:`peregrine.iqgen.bits.amplitude_factory` module contains classes and +functions related to object factory for amplitude objects. + +""" + +from peregrine.iqgen.bits.amplitude_poly import AmplitudePoly as PolyAmplitude +from peregrine.iqgen.bits.amplitude_sine import AmplitudeSine as SineAmplitude + + +class ObjectFactory(object): + ''' + Object factory for amplitude objects. + ''' + + def __init__(self): + super(ObjectFactory, self).__init__() + + def toMapForm(self, obj): + t = type(obj) + if t is PolyAmplitude: + return self.__PolyAmplitude_ToMap(obj) + elif t is SineAmplitude: + return self.__SineAmplitude_ToMap(obj) + else: + raise ValueError("Invalid object type") + + def fromMapForm(self, data): + t = data['type'] + if t == 'PolyAmplitude': + return self.__MapTo_PolyAmplitude(data) + elif t == 'SineAmplitude': + return self.__MapTo_SineAmplitude(data) + else: + raise ValueError("Invalid object type") + + def __PolyAmplitude_ToMap(self, obj): + data = {'type': 'PolyAmplitude', 'coeffs': obj.coeffs} + return data + + def __SineAmplitude_ToMap(self, obj): + data = {'type': 'SineAmplitude', + 'initial': obj.initial, + 'amplitude': obj.amplitude, + 'period': obj.period_s} + return data + + def __MapTo_PolyAmplitude(self, data): + coeffs = data['coeffs'] + return PolyAmplitude(coeffs) + + def __MapTo_SineAmplitude(self, data): + initial = data['initial'] + amplitude = data['amplitude'] + period = data['period'] + return SineAmplitude(initial, amplitude, period) + +factoryObject = ObjectFactory() diff --git a/peregrine/iqgen/bits/amplitude_poly.py b/peregrine/iqgen/bits/amplitude_poly.py new file mode 100644 index 0000000..0f9275a --- /dev/null +++ b/peregrine/iqgen/bits/amplitude_poly.py @@ -0,0 +1,95 @@ +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Valeri Atamaniouk +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + +""" +The :mod:`peregrine.iqgen.amplitude_poly` module contains classes and functions +related to implementation of polynomial-based amplitude class. + +""" + +from peregrine.iqgen.bits.amplitude_base import AmplitudeBase + +import numpy + + +class AmplitudePoly(AmplitudeBase): + ''' + Amplitude control with polynomial dependency over time. + ''' + + def __init__(self, coeffs): + ''' + Constructs polynomial amplitude control object. + + Parameters + coeffs : array-like + Polynomial coefficients + ''' + super(AmplitudePoly, self).__init__() + + self.coeffs = tuple([x for x in coeffs]) + if len(coeffs) > 0: + self.poly = numpy.poly1d(coeffs) + else: + self.poly = None + + def __str__(self): + ''' + Constructs literal presentation of object. + + Returns + ------- + string + Literal presentation of object + ''' + return "AmplitudePoly(c={})".format(self.coeffs) + + def applyAmplitude(self, signal, userTimeAll_s): + ''' + Applies amplitude modulation to signal. + + This method applies polynomial modulation. + + Parameters + ---------- + signal : numpy.ndarray + Signal sample vector. Each element defines signal amplitude in range + [-1; +1]. This vector is modified in place. + userTimeAll_s : numpy.ndarray + Sample time vector. Each element defines sample time in seconds. + + Returns + ------- + numpy.ndarray + Array with output samples + ''' + + poly = self.poly + if poly is not None: + amplitudeVector = poly(userTimeAll_s) + signal *= amplitudeVector + + return signal + + def computeMeanPower(self): + ''' + Computes mean signal power. + + Returns + ------- + float + Mean signal power for the configured amplitude + ''' + poly = self.poly + if poly is not None: + result = numpy.square(poly(0.)) + else: + result = 1. + return result diff --git a/peregrine/iqgen/bits/amplitude_sine.py b/peregrine/iqgen/bits/amplitude_sine.py new file mode 100644 index 0000000..ee18452 --- /dev/null +++ b/peregrine/iqgen/bits/amplitude_sine.py @@ -0,0 +1,91 @@ +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Valeri Atamaniouk +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + +""" +The :mod:`peregrine.iqgen.amplitude_sine` module contains classes and functions +related to implementation of sine-based amplitude class. + +""" + +from peregrine.iqgen.bits.amplitude_base import AmplitudeBase + +import numpy +import scipy.constants + + +class AmplitudeSine(AmplitudeBase): + ''' + Amplitude control with sine modulation over time. + ''' + + def __init__(self, initial, amplitude, period_s): + ''' + Constructs sine amplitude control object. + + Parameters + initial : float + Initial amplitude value (median) + amplitude : float + Amplitude of change + period_s : float + Period of change in seconds + ''' + super(AmplitudeSine, self).__init__() + self.initial = initial + self.amplitude = amplitude + self.period_s = period_s + self.c = 2. * scipy.constants.pi / period_s + + def __str__(self): + ''' + Constructs literal presentation of object. + + Returns + ------- + string + Literal presentation of object + ''' + return "AmplitudeSine(base={}, amp={}, p={} s)".\ + format(self.initial, self.amplitude, self.period_s) + + def applyAmplitude(self, signal, userTimeAll_s): + ''' + Applies amplitude modulation to signal. + + Parameters + ---------- + signal : numpy.ndarray + Signal sample vector. Each element defines signal amplitude in range + [-1; +1]. This vector is modified in place. + userTimeAll_s : numpy.ndarray + Sample time vector. Each element defines sample time in seconds. + + Returns + ------- + numpy.ndarray + Array with output samples + ''' + + ampAll = numpy.sin(userTimeAll_s * self.c) * self.amplitude + self.initial + + return numpy.multiply(signal, ampAll, out=signal) + + def computeMeanPower(self): + ''' + Computes mean signal power. + + Returns + ------- + float + Mean signal power for the configured amplitude + ''' + amplitude = self.amplitude * 0.707 + self.initial + result = numpy.square(amplitude) + return result diff --git a/peregrine/iqgen/bits/doppler_base.py b/peregrine/iqgen/bits/doppler_base.py new file mode 100644 index 0000000..ebddc8e --- /dev/null +++ b/peregrine/iqgen/bits/doppler_base.py @@ -0,0 +1,334 @@ +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Valeri Atamaniouk +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + +""" +The :mod:`peregrine.iqgen.doppler_base` module contains classes and functions +related to base implementation of doppler class. + +""" + +import scipy.constants +import numpy + + +class DopplerBase(object): + ''' + Doppler control for a signal source that moves with a constant speed. + ''' + + def __init__(self, distance0_m=0., tec_epm2=50., dtype=numpy.longdouble): + ''' + Constructs doppler base object for movement control. + + Parameters + ---------- + distance0_m : float + Distance to object in meters at time 0. + tec_epm2 : float + Total free electron content integrated along line of sight to the object + in electrons per m^2. + dtype : object, optional + Numpy type for sample computations. + ''' + super(DopplerBase, self).__init__() + self.distance0_m = distance0_m + self.tec_epm2 = tec_epm2 + self.dtype = dtype + self.codeDopplerIgnored = False + self.twoPi = scipy.constants.pi * 2. + + def isCodeDopplerIgnored(self): + ''' + Checks if the code is ignoring doppler control. + + Returns + ------- + bool + When True, the sample generator ignores doppler shift for data and code + processing. + ''' + return self.codeDopplerIgnored + + def setCodeDopplerIgnored(self, value): + ''' + Changes doppler control for data and code processing + + Parameters + ---------- + value : bool + True - ignore doppler for code and data processing, False - apply doppler. + ''' + self.codeDopplerIgnored = value + + def computeSignalDelayS(self, frequency_hz): + ''' + Computes delay in seconds for an epoch time (time 0) for a given carrier + frequency. + + The method computes signal delay, which is a sum of the following + parameters: + - Distance to object divided per speed of light + - Ionospheric delay according to TEC value for the given frequency + - Tropospheric delay + + Parameters + ---------- + frequency_hz : float + Signal frequency in hertz. + + Returns + ------- + float + Signal delay in seconds. + ''' + distanceDelay_s = self.distance0_m / scipy.constants.c + ionoDelay_s = 40.3 * self.tec_epm2 / numpy.square(frequency_hz) + delay_s = distanceDelay_s + ionoDelay_s + return delay_s + + def applySignalDelays(self, userTimeAll_s, carrierSignal): + ''' + Modifies time vector in accordance to signal delays due to distance and + atmospheric delays. + + Parameters + ---------- + userTimeAll_s : numpy.ndvector(shape=(n), dtype=numpy.float) + Vector of time stamps for which samples are generated + carrierSignal : object + Signal parameters object + + Returns + ------- + numpy.ndvector(shape=(n), dtype=numpy.float) + Vector of sample time stamps updated according to computed delays. + ''' + signalDelay_s = self.computeSignalDelayS(carrierSignal.CENTER_FREQUENCY_HZ) + return userTimeAll_s - signalDelay_s + + def computeDistanceM(self, svTime_s): + ''' + Computes doppler shift in meters. + + Parameters + ---------- + svTime_s : float + Time in seconds at which distance is computed. Please note that is not + a time of the observer. + + Returns + ------- + float + Distance to satellite in meters. + ''' + raise NotImplementedError() + + def computeSpeedMps(self, svTime_s): + ''' + Computes speed along the vecto2r to satellite in meters per second. + + Parameters + ---------- + svTime_s : float + Time in seconds at which speed is computed. Please note that is not + a time of the observer. + + Returns + ------- + float + Speed of satellite in meters per second. + ''' + raise NotImplementedError() + + def computeBatch(self, + userTimeAll_s, + amplitude, + carrierSignal, + ifFrequency_hz, + message, + code, + outputConfig, + debug): + ''' + Computes signal samples for the doppler object. + + Parameters + ---------- + userTimeAll_s : numpy.ndarray(dtype=numpy.float) + Sample timestamps in seconds + amplitude : float + Signal amplitude object. + carrierSignal : object + Carrier frequency object + ifFrequency_hz: float + Intermediate frequency in hertz + message : object + Message object for providing access to symbols + code : object + PRN code object for providing access to chips + debug : bool + Debug flag + + Returns + ------- + signal : numpy.ndarray(n_samples, dtype=float) + Generated samples + userTimeX_s : float + End of interval time in seconds + chipAll_idx : numpy.ndarray(n_samples, dtype=float) + Code chip phases for the samples + chips : numpy.ndarray(n_samples, dtype=int) + Code combined with data + ''' + + userTimeAll_s = self.applySignalDelays(userTimeAll_s, carrierSignal) + + # Computing doppler coefficients + twoPi = self.twoPi + + # Sine wave phase without doppler + phaseAll = userTimeAll_s * (twoPi * ifFrequency_hz) + # Get doppler shift in meters + doppler_m = self.computeDopplerShiftM(userTimeAll_s) + # Doppler for carrier center frequency + carrFreqRatio = -carrierSignal.CENTER_FREQUENCY_HZ / scipy.constants.c + phaseAll += doppler_m * (carrFreqRatio * twoPi) + + # Convert phase to signal value and multiply by amplitude + signal = scipy.cos(phaseAll) + + if amplitude: + amplitude.applyAmplitude(signal, userTimeAll_s) + + # PRN and data index computation + chipAll_idx = userTimeAll_s * carrierSignal.CODE_CHIP_RATE_HZ + if self.codeDopplerIgnored: + pass + else: + # Computing doppler coefficients + chipFreqRatio = -carrierSignal.CODE_CHIP_RATE_HZ / scipy.constants.c + chipAll_idx += doppler_m * chipFreqRatio + + chips = self.computeDataNChipVector(chipAll_idx, + carrierSignal, + message, + code) + + # Combine data and sine wave + signal *= chips + + # Generate debug data + doppler_hz = self.computeDopplerShiftHz(userTimeAll_s, + carrierSignal) if debug else None + return (signal, doppler_hz) + + @staticmethod + def computeDeltaUserTimeS(userTime0_s, n_samples, outputConfig): + ''' + Helper for computing generation interval duration in seconds. + + Parameters + ---------- + userTime0_s : float + Generation interval start + n_samples : int + Number of samples in the generation interval + outputConfig : object + Output configuration. + + Returns + ------- + float + Generation interval duration in seconds + ''' + deltaUserTime_s = float(n_samples) / outputConfig.SAMPLE_RATE_HZ + return deltaUserTime_s + + @staticmethod + def computeDopplerHz(frequency_hz, speed_mps): + ''' + Generic method for doppler shift computation. + + Parameters + ---------- + frequency_hz : float + Frequency in hertz for which doppler is computed. + speed_mps : float + Speed in meters per second for which doppler is computed. + + Returns + ------- + float + Doppler shift value in hertz. + ''' + doppler_hz = -frequency_hz / scipy.constants.c * speed_mps + return doppler_hz + + def computeDataNChipVector(self, chipAll_idx, carrierSignal, message, code): + ''' + Helper for computing vector that combines data and code chips. + + Parameters + ---------- + chipAll_idx : ndarray + vector of chip phases + carrierSignal : object + Signal description object + messge : object + Data bits source + code : objects + Code chips source + + Returns + ------- + ndarray + Array of code chips multiplied with data bits + ''' + + chipAll_long = chipAll_idx.astype(numpy.long) + dataBits = message.getDataBits( + chipAll_long / carrierSignal.CHIP_TO_SYMBOL_DIVIDER) + result = code.combineData(chipAll_long, dataBits) + + return result + + def computeDopplerShiftM(self, userTimeAll_s): + ''' + Method to compute metric doppler shift + + Parameters + ---------- + userTimeAll_s : numpy.ndarray(shape=(1, nSamples), dtype=numpy.float) + Time vector for sample timestamps in seconds + + Returns + ------- + numpy.ndarray(shape=(1, nSamples), dtype=numpy.float) + Computed doppler shift in meters + ''' + raise NotImplementedError("Metric doppler computation is not implemented") + + def computeDopplerShiftHz(self, userTimeAll_s, carrierSignal): + ''' + Method to compute doppler shift in Hz. + + Parameters + ---------- + userTimeAll_s : numpy.ndarray(shape=(1, nSamples), dtype=numpy.float) + Time vector for sample timestamps in seconds + carrierSignal : object + Carrier signal parameters + + Returns + ------- + numpy.ndarray(shape=(1, nSamples), dtype=numpy.float) + Computed doppler frquency shift in hertz + ''' + raise NotImplementedError("Metric doppler computation is not implemented") diff --git a/peregrine/iqgen/bits/doppler_factory.py b/peregrine/iqgen/bits/doppler_factory.py new file mode 100644 index 0000000..6de2a35 --- /dev/null +++ b/peregrine/iqgen/bits/doppler_factory.py @@ -0,0 +1,83 @@ +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Valeri Atamaniouk +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + +""" +The :mod:`peregrine.iqgen.bits.doppler_factory` module contains classes and +functions related to object factory for doppler control objects. + +""" + +from peregrine.iqgen.bits.doppler_poly import Doppler as PolyDoppler +from peregrine.iqgen.bits.doppler_sine import Doppler as SineDoppler + + +class ObjectFactory(object): + ''' + Object factory for doppler control objects. + ''' + + def __init__(self): + super(ObjectFactory, self).__init__() + + def toMapForm(self, obj): + t = type(obj) + if t is PolyDoppler: + return self.__PolyDoppler_ToMap(obj) + elif t is SineDoppler: + return self.__SineDoppler_ToMap(obj) + else: + raise ValueError("Invalid object type") + + def fromMapForm(self, data): + t = data['type'] + if t == 'PolyDoppler': + return self.__MapTo_PolyDoppler(data) + elif t == 'SineDoppler': + return self.__MapTo_SineDoppler(data) + else: + raise ValueError("Invalid object type") + + def __PolyDoppler_ToMap(self, obj): + data = {'type': 'PolyDoppler', + 'distance0_m': obj.distance0_m, + 'tec_epm2': obj.tec_epm2, + 'coeffs': obj.coeffs} + return data + + def __SineDoppler_ToMap(self, obj): + data = {'type': 'SineDoppler', + 'distance0_m': obj.distance0_m, + 'tec_epm2': obj.tec_epm2, + 'speed0_mps': obj.speed0_mps, + 'amplutude_mps': obj.amplutude_mps, + 'period_s': obj.period_s} + return data + + def __MapTo_PolyDoppler(self, data): + distance0_m = data['distance0_m'] + tec_epm2 = data['tec_epm2'] + coeffs = data['coeffs'] + return PolyDoppler(distance0_m=distance0_m, + tec_epm2=tec_epm2, + coeffs=coeffs) + + def __MapTo_SineDoppler(self, data): + distance0_m = data['distance0_m'] + tec_epm2 = data['tec_epm2'] + speed0_mps = data['speed0_mps'] + amplutude_mps = data['amplutude_mps'] + period_s = data['period_s'] + return SineDoppler(distance0_m=distance0_m, + tec_epm2=tec_epm2, + speed0_mps=speed0_mps, + amplutude_mps=amplutude_mps, + period_s=period_s) + +factoryObject = ObjectFactory() diff --git a/peregrine/iqgen/bits/doppler_poly.py b/peregrine/iqgen/bits/doppler_poly.py new file mode 100644 index 0000000..d897716 --- /dev/null +++ b/peregrine/iqgen/bits/doppler_poly.py @@ -0,0 +1,246 @@ +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Valeri Atamaniouk +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + +""" +The :mod:`peregrine.iqgen.bits.oppler_poly` module contains classes and +functions related to generation of signals with polynomial-based movement. + +""" + +import numpy +import scipy.constants +from peregrine.iqgen.bits.doppler_base import DopplerBase + + +class Doppler(DopplerBase): + ''' + Doppler control for an object that has constant acceleration. Such signal has + constant doppler value with a possible sign invert. + ''' + + TWO_PI = scipy.constants.pi * 2 + + def __init__(self, distance0_m, tec_epm2, coeffs): + ''' + Constructs doppler control object for linear acceleration. + + Parameters + ---------- + distance0_m : float + Distance to object in meters at time 0. + tec_epm2 : float + Total free electron content integrated along line of sight to the object + in electrons per m^2. + coeffs : array-like + Phase shift coefficients. Phase chift will be computed as: + C_n*t^n + C_(n-1)^(n-1) + ... + C_2*t^2 + C_1*t + C_0 + C_n..C_0 - values for speed of light + ''' + super(Doppler, self).__init__(distance0_m=distance0_m, + tec_epm2=tec_epm2) + self.coeffs = tuple([x for x in coeffs]) + self.n_coeffs = len(coeffs) + self.speedPoly = None + self.distancePoly = None + if self.n_coeffs > 0: + new_coeffs = [] + self.n_coeffs += 1 + for idx, c in enumerate(coeffs): + order = self.n_coeffs - idx - 1 + new_coeffs.append(c / order) + new_coeffs.append(0.) + self.distancePoly = numpy.poly1d(new_coeffs) + self.distanceCoeffs = new_coeffs + if self.n_coeffs > 1: + self.speedPoly = numpy.poly1d(coeffs) + else: + self.distanceCoeffs = None + + def __str__(self): + ''' + Constructs literal presentation of object. + + Returns + ------- + string + Literal presentation of object + ''' + return "DopplerPoly(coeffs={}, distance0_m={}," \ + " tec_epm2={} codeDopplerIgnored={})". \ + format(self.coeffs, self.distance0_m, + self.tec_epm2, self.codeDopplerIgnored) + + def computeDistanceM(self, svTime_s): + ''' + Computes doppler shift in meters. + + Parameters + ---------- + svTime_s : float + Time in seconds at which distance is computed. Please note that is not + a time of the observer. + + Returns + ------- + float + Distance to satellite in meters. + ''' + poly = self.distancePoly + if poly is not None: + return poly(svTime_s) # self.coeffs[cnt - 1] + else: + return 0. + + def computeSpeedMps(self, svTime_s): + ''' + Computes speed along the vector to satellite in meters per second. + + Parameters + ---------- + svTime_s : float + Time in seconds at which speed is computed. Please note that is not + a time of the observer. + + Returns + ------- + float + Speed of satellite in meters per second. + ''' + poly = self.speedPoly + if poly is not None: + return poly(svTime_s) + else: + return 0. + + def computeDopplerShiftM(self, userTimeAll_s): + ''' + Method to compute metric doppler shift + + Parameters + ---------- + userTimeAll_s : numpy.ndarray(shape=(1, nSamples), dtype=numpy.float) + Time vector for sample timestamps in seconds + + Returns + ------- + numpy.ndarray(shape=(1, nSamples), dtype=numpy.float) + Computed doppler shift in meters + ''' + distancePoly = self.distancePoly + if distancePoly is not None: + # Slower, but simple + doppler_m = distancePoly(userTimeAll_s) + else: + # No phase shift + doppler_m = numpy.zeros_like(userTimeAll_s) + return doppler_m + + def computeDopplerShiftHz(self, userTimeAll_s, carrierSignal): + ''' + Method to compute doppler shift in Hz. + + Parameters + ---------- + userTimeAll_s : numpy.ndarray(shape=(1, nSamples), dtype=numpy.float) + Time vector for sample timestamps in seconds + carrierSignal : object + Carrier signal parameters + + Returns + ------- + numpy.ndarray(shape=(1, nSamples), dtype=numpy.float) + Computed doppler frquency shift in hertz + ''' + speedPoly = self.speedPoly + if speedPoly is not None: + # Slower, but simple + c0 = -carrierSignal.CENTER_FREQUENCY_HZ / scipy.constants.c + doppler_hz = speedPoly(userTimeAll_s) * c0 + else: + # No phase shift + doppler_hz = numpy.zeros_like(userTimeAll_s) + return doppler_hz + + +def linearDoppler(distance0_m, + tec_epm2, + frequency_hz, + doppler0_hz, + dopplerChange_hzps): + ''' + Makes an object that corresponds to linear doppler change. + + Parameters + ---------- + distance0_m : float + Initial distance to object. + doppler0_hz : float + Initial doppler shift in hz. + frequency_hz + Carrier frequency in Hz. + dopplerChange_hzps : float + Doppler shift rate in Hz per second. + + Returns + ------- + Doppler + object that implments constant acceleration logic. + ''' + speed0_mps = -scipy.constants.c / frequency_hz * doppler0_hz + accel_mps2 = -scipy.constants.c / frequency_hz * dopplerChange_hzps + + return Doppler(distance0_m=distance0_m, + tec_epm2=tec_epm2, + coeffs=(accel_mps2, speed0_mps)) + + +def constDoppler(distance0_m, tec_epm2, frequency_hz, doppler_hz): + ''' + Makes an object that corresponds to a constant doppler value. + + Parameters + ---------- + distance0_m : float + Initial distance to object. + frequency_hz : float + Carrier frequency in Hz. + doppler_hz : float + Doppler shift in Hz. + + Returns + ------- + Doppler + Object that implements constant speed logic. + ''' + speed_mps = -scipy.constants.c / frequency_hz * doppler_hz + return Doppler(distance0_m=distance0_m, tec_epm2=tec_epm2, coeffs=(speed_mps,)) + + +def zeroDoppler(distance_m, tec_epm2, frequency_hz): + ''' + Makes an object that corresponds to zero doppler change. + + Parameters + ---------- + distance0_m : float + Initial distance to object. + doppler0_hz : float + Initial doppler shift in hz. + frequency_hz + Carrier frequency in Hz. + dopplerChange_hzps : float + Doppler shift rate in Hz per second. + + Returns + ------- + Doppler + object that implments constant acceleration logic. + ''' + return Doppler(distance0_m=distance_m, tec_epm2=tec_epm2, coeffs=()) diff --git a/peregrine/iqgen/bits/doppler_sine.py b/peregrine/iqgen/bits/doppler_sine.py new file mode 100644 index 0000000..95d85e3 --- /dev/null +++ b/peregrine/iqgen/bits/doppler_sine.py @@ -0,0 +1,199 @@ +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Valeri Atamaniouk +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. +""" +The :mod:`peregrine.iqgen.bits.doppler_sine` module contains classes and +functions related to generation of signals with circular changing doppler. + +""" + +from peregrine.iqgen.bits.doppler_base import DopplerBase + +import scipy.constants +import numpy + + +class Doppler(DopplerBase): + ''' + Doppler control for an object that has peridic acceleration. + ''' + + TWO_PI = scipy.constants.pi * 2. + + def __init__(self, distance0_m, tec_epm2, speed0_mps, amplutude_mps, period_s): + ''' + Constructs doppler control object for linear acceleration. + + Parameters + ---------- + distance0_m : float + Distance to object in meters at time 0. + tec_epm2 : float + Total free electron content integrated along line of sight to the object + in electrons per m^2. + speed0_mps : float + Speed of satellite at time 0 in meters per second. + amplutude_mps : float + Amplitude of change + period_s : float + Period of change + ''' + super(Doppler, self).__init__(distance0_m=distance0_m, + tec_epm2=tec_epm2) + self.speed0_mps = speed0_mps + self.amplutude_mps = amplutude_mps + self.period_s = period_s + + def __str__(self): + ''' + Constructs literal presentation of object. + + Returns + ------- + string + Literal presentation of object + ''' + return "SineDoppler(distance0_m={}, tec_epm2={}," \ + " speed0_mps={}, amplitude_mps={}, period_s={}," \ + " codeDopplerIgnored={})".\ + format(self.distance0_m, self.tec_epm2, self.speed0_mps, + self.amplutude_mps, self.period_s, self.codeDopplerIgnored) + + def __repr__(self): + ''' + Constructs python expression presentation of object. + + Returns + ------- + string + Python expression presentation of object + ''' + return "Doppler({}, {}, {}, {}, {})".format(self.distance0_m, + self.tec_epm2, + self.speed0_mps, + self.amplutude_mps, + self.period_s) + + def computeDistanceM(self, svTime_s): + ''' + Computes doppler shift in meters. + + Parameters + ---------- + svTime_s : float + Time in seconds at which distance is computed. Please note that is not + a time of the observer. + + Returns + ------- + float + Distance to satellite in meters. + ''' + return self.distance0_m + self.speed0_mps * svTime_s + \ + self.amplutude_mps * \ + (1 - numpy.cos(Doppler.TWO_PI * svTime_s / self.period_s)) + + def computeSpeedMps(self, svTime_s): + ''' + Computes speed along the vector to satellite in meters per second. + + Parameters + ---------- + svTime_s : float + Time in seconds at which speed is computed. Please note that is not + a time of the observer. + + Returns + ------- + float + Speed of satellite in meters per second. + ''' + return self.speed0_mps + self.amplutude_mps * \ + numpy.sin(Doppler.TWO_PI * svTime_s / self.period_s) + + def computeDopplerShiftM(self, userTimeAll_s): + ''' + Method to compute metric doppler shift + + Parameters + ---------- + userTimeAll_s : numpy.ndarray(shape=(1, nSamples), dtype=numpy.float) + Time vector for sample timestamps in seconds + + Returns + ------- + numpy.ndarray(shape=(1, nSamples), dtype=numpy.float) + Computed doppler shift in meters + ''' + D_0 = self.speed0_mps + D_1 = self.amplutude_mps * self.period_s / self.twoPi + D_2 = self.twoPi / self.period_s + + doppler_m = numpy.cos(D_2 * userTimeAll_s) + doppler_m -= 1. + doppler_m *= -D_1 + if D_0: + doppler_m += D_0 * userTimeAll_s + + return doppler_m + + def computeDopplerShiftHz(self, userTimeAll_s, carrierSignal): + ''' + Method to compute doppler shift in Hz. + + Parameters + ---------- + userTimeAll_s : numpy.ndarray(shape=(1, nSamples), dtype=numpy.float) + Time vector for sample timestamps in seconds + carrierSignal : object + Carrier signal parameters + + Returns + ------- + numpy.ndarray(shape=(1, nSamples), dtype=numpy.float) + Computed doppler frquency shift in hertz + ''' + D_0 = self.speed0_mps + D_1 = self.amplutude_mps + D_2 = self.twoPi / self.period_s + + doppler_hz = numpy.sin(D_2 * userTimeAll_s) * D_1 + if D_0: + doppler_hz += D_0 + doppler_hz *= -carrierSignal.CENTER_FREQUENCY_HZ / scipy.constants.c + return doppler_hz + + +def sineDoppler(distance0_m, tec_epm2, frequency_hz, doppler0_hz, dopplerAmplitude_hz, dopplerPeriod_s): + ''' + Makes an object that corresponds to linear doppler change. + + Parameters + ---------- + distance0_m : float + Initial distance to object. + frequency_hz + Carrier frequency in Hz. + doppler0_hz : float + Initial doppler shift in hz. + dopplerAmplitude_hz : float + Doppler change amplitude in Hz + dopplerPeriod_s : float + Doppler change period in seconds + + Returns + ------- + Doppler + object that implments constant acceleration logic. + ''' + dopplerCoeff = -scipy.constants.c / frequency_hz + speed0_mps = dopplerCoeff * doppler0_hz + amplitude_mps = dopplerCoeff * dopplerAmplitude_hz + + return Doppler(distance0_m, tec_epm2, speed0_mps, amplitude_mps, dopplerPeriod_s) diff --git a/peregrine/iqgen/bits/encoder_1bit.py b/peregrine/iqgen/bits/encoder_1bit.py new file mode 100644 index 0000000..2690bcc --- /dev/null +++ b/peregrine/iqgen/bits/encoder_1bit.py @@ -0,0 +1,82 @@ +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Valeri Atamaniouk +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + +""" +The :mod:`peregrine.iqgen.bits.encoder_1bit` module contains classes and +functions related to generating single bit signal output. + +""" + +from peregrine.iqgen.bits.encoder_base import Encoder + + +class BandBitEncoder(Encoder): + ''' + Base class for single bit encoding. + ''' + + def __init__(self, bandIndex): + ''' + Initializes encoder object. + + Parameters + ---------- + bandIndex : int + Index of the band in the generated sample matrix. + ''' + super(BandBitEncoder, self).__init__() + self.bandIndex = bandIndex + + def addSamples(self, sample_array): + ''' + Extracts samples of the supported band and coverts them into bit stream. + + Parameters + ---------- + sample_array : numpy.ndarray((4, N)) + Sample vectors ordered by band index. + + Returns + ------- + numpy.ndarray((N), dtype=numpy.uint8) + Array of type uint8 containing the encoded data. + ''' + band_samples = sample_array[self.bandIndex] + n_samples = len(band_samples) + + self.ensureExtraCapacity(n_samples) + start = self.n_bits + end = start + n_samples + self.bits[start:end] = BandBitEncoder.convertBand(band_samples) + self.n_bits = end + + if (self.n_bits >= Encoder.BLOCK_SIZE): + return self.encodeValues() + + return Encoder.EMPTY_RESULT + + @staticmethod + def convertBand(band_samples): + ''' + Helper method for converting sampled signal band into output bits. + + The samples are compared to 0. Positive values yield value of False. + + Parameters + ---------- + band_samples : numpy.ndarray((N)) + Vector of signal samples + + Returns + ------- + signs : numpy.ndarray((N), dtype=numpy.bool) + Boolean vector of sample signs + ''' + return band_samples < 0 diff --git a/peregrine/iqgen/bits/encoder_2bits.py b/peregrine/iqgen/bits/encoder_2bits.py new file mode 100644 index 0000000..286a884 --- /dev/null +++ b/peregrine/iqgen/bits/encoder_2bits.py @@ -0,0 +1,117 @@ +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Valeri Atamaniouk +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + +""" +The :mod:`peregrine.iqgen.bits.encoder_2bits` module contains classes and +functions related to generating two bits signal output. + +""" + +import numpy + +from peregrine.iqgen.bits.encoder_base import Encoder + + +class BandTwoBitsEncoder(Encoder): + ''' + Base class for two bits encoding. + ''' + + def __init__(self, bandIndex): + ''' + Initializes encoder object. + + Parameters + ---------- + bandIndex : int + Index of the band in the generated sample matrix. + ''' + super(BandTwoBitsEncoder, self).__init__() + self.bandIndex = bandIndex + + @staticmethod + def convertBand(band_samples): + ''' + Helper method for converting sampled signal band into output bits. + + For the sign, the samples are compared to 0. Positive values yield sign of + True. + + The method builds a power histogram from signal samples. After a histogram + is built, the 67% power boundary is located. All samples, whose power is + lower, than the boundary, are reported as False. + + Parameters + ---------- + band_samples : numpy.ndarray + Vector of signal samples + + Returns + ------- + signs : numpy.ndarray(dtype=numpy.bool) + Boolean vector of sample signs: True for positive, False for negative + amps : numpy.ndarray(dtype=numpy.bool) + Boolean vector of sample power: True for high power, False for low power + ''' + + # Signal power is a square of the amplitude + power = numpy.square(band_samples) + totalPower = numpy.sum(power) + totalPowerLimit = totalPower * 0.67 + + # Build histogram to find 67% power + totalBins = 30 + hist, edges = numpy.histogram(power, + bins=totalBins, + density=False) + avg = (edges[:-1] + edges[1:]) * 0.5 + powers = numpy.cumsum(hist * avg) + idx = numpy.searchsorted(powers, totalPowerLimit, side="right") + powerLimit = avg[idx] + + # Signal sign + signs = band_samples > 0 + amps = power >= powerLimit + + return signs, amps + + def addSamples(self, sample_array): + ''' + Extracts samples of the supported band and coverts them into bit stream. + + Parameters + ---------- + sample_array : numpy.ndarray((4, N)) + Sample vectors ordered by band index. + + Returns + ------- + numpy.ndarray(dtype=numpy.uint8) + Array of type uint8 containing the encoded data. + ''' + band_samples = sample_array[self.bandIndex] + n_samples = len(band_samples) + + # Signal signs and amplitude + signs, amps = self.convertBand(band_samples) + + self.ensureExtraCapacity(n_samples * 2) + + bits = self.bits + start = self.n_bits + end = start + n_samples * 2 + bits[start + 0:end:2] = signs + bits[start + 1:end:2] = amps + self.n_bits = end + + if (self.n_bits >= Encoder.BLOCK_SIZE): + return self.encodeValues() + else: + return Encoder.EMPTY_RESULT diff --git a/peregrine/iqgen/bits/encoder_base.py b/peregrine/iqgen/bits/encoder_base.py new file mode 100644 index 0000000..ee6b134 --- /dev/null +++ b/peregrine/iqgen/bits/encoder_base.py @@ -0,0 +1,105 @@ +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Valeri Atamaniouk +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + +""" +The :mod:`peregrine.iqgen.bits.encoder_base` module contains classes and +functions related to generating signal output. + +""" + +import numpy + + +class Encoder(object): + ''' + Base encode class. + + Encoder accepts sequence of signal arrays as input and produces byte arrays + as output. + ''' + # Block size for encode. Must be a multiple of 8. + BLOCK_SIZE = 1024 * 8 + + EMPTY_RESULT = numpy.ndarray(0, dtype=numpy.uint8) # Internal empty array + + def __init__(self, bufferSize=1000): + ''' + Constructs encoder. + + Parameters + ---------- + bufferSize : int, optional + Size of the internal buffer to batch-process samples + ''' + self.bits = numpy.ndarray(bufferSize, dtype=numpy.int8) + self.n_bits = 0 + + def addSamples(self, sample_array): + ''' + Extracts samples of the supported band and coverts them into bit stream. + + Parameters + ---------- + sample_array : numpy.ndarray((4, N)) + Sample vectors ordered by band index. + + Returns + ------- + numpy.ndarray(dtype=numpy.uint8) + Array of type uint8 containing the encoded data. + ''' + return Encoder.EMPTY_RESULT + + def flush(self): + ''' + Flushes the data in the buffer. + + Returns + ------- + ndarray + Array of type uint8 containing the encoded data. + ''' + if self.n_bits > 0 and self.n_bits % 8 != 0: + self.bits += (8 - self.n_bits % 8) + res = numpy.packbits(self.bits[0:self.n_bits]) + self.n_bits = 0 + return res + + def encodeValues(self): + ''' + Converts buffered bit data into packed array. + + The method coverts multiple of 8 bits into single output byte. + + Returns + ------- + ndarray + Array of type uint8 containing the encoded data. + ''' + n_bytes = self.n_bits / 8 + n_offset = n_bytes * 8 + n_left = self.n_bits - n_offset + res = numpy.packbits(self.bits[0: n_offset]) + self.bits[0:n_left] = self.bits[n_offset:n_offset + n_left] + self.n_bits = n_left + return res + + def ensureExtraCapacity(self, extraBits): + ''' + Method verifies that current array has sufficient capacity to hold + additional bits. + + Parameters + ---------- + extraBits : int + Number of extra bits to reserve space for + ''' + if len(self.bits) < self.n_bits + extraBits: + self.bits.resize(self.n_bits + extraBits) diff --git a/peregrine/iqgen/bits/encoder_factory.py b/peregrine/iqgen/bits/encoder_factory.py new file mode 100644 index 0000000..ab73321 --- /dev/null +++ b/peregrine/iqgen/bits/encoder_factory.py @@ -0,0 +1,109 @@ +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Valeri Atamaniouk +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + +""" +The :mod:`peregrine.iqgen.bits.encoder_factory` module contains classes and +functions related to object factory for output encoder objects. + +""" + +from peregrine.iqgen.bits.encoder_gps import GPSL1BitEncoder +from peregrine.iqgen.bits.encoder_gps import GPSL2BitEncoder +from peregrine.iqgen.bits.encoder_gps import GPSL1L2BitEncoder +from peregrine.iqgen.bits.encoder_gps import GPSL1TwoBitsEncoder +from peregrine.iqgen.bits.encoder_gps import GPSL2TwoBitsEncoder +from peregrine.iqgen.bits.encoder_gps import GPSL1L2TwoBitsEncoder + + +class ObjectFactory(object): + ''' + Object factory for encoder objects. + ''' + + def __init__(self): + super(ObjectFactory, self).__init__() + + def toMapForm(self, obj): + t = type(obj) + if t is GPSL1BitEncoder: + return self.__GPSL1BitEncoder_ToMap(obj) + elif t is GPSL2BitEncoder: + return self.__GPSL2BitEncoder_ToMap(obj) + elif t is GPSL1L2BitEncoder: + return self.__GPSL1L2BitEncoder_ToMap(obj) + elif t is GPSL1TwoBitsEncoder: + return self.__GPSL1TwoBitsEncoder_ToMap(obj) + elif t is GPSL2TwoBitsEncoder: + return self.__GPSL2TwoBitsEncoder_ToMap(obj) + elif t is GPSL1L2TwoBitsEncoder: + return self.__GPSL1L2TwoBitsEncoder_ToMap(obj) + else: + raise ValueError("Invalid object type") + + def fromMapForm(self, data): + t = data['type'] + if t == 'GPSL1BitEncoder': + return self.__MapTo_GPSL1BitEncoder(data) + elif t == 'GPSL2BitEncoder': + return self.__MapTo_GPSL2BitEncoder(data) + elif t == 'GPSL2BitEncoder': + return self.__MapTo_GPSL1L2BitEncoder(data) + elif t == 'GPSL1TwoBitsEncoder': + return self.__MapTo_GPSL1TwoBitsEncoder(data) + elif t == 'GPSL2TwoBitsEncoder': + return self.__MapTo_GPSL2TwoBitsEncoder(data) + elif t == 'GPSL1L2TwoBitsEncoder': + return self.__MapTo_GPSL1L2TwoBitsEncoder(data) + else: + raise ValueError("Invalid object type") + + def __GPSL1BitEncoder_ToMap(self, obj): + data = {'type': 'GPSL1BitEncoder'} + return data + + def __GPSL2BitEncoder_ToMap(self, obj): + data = {'type': 'GPSL2BitEncoder'} + return data + + def __GPSL1L2BitEncoder_ToMap(self, obj): + data = {'type': 'GPSL1L2BitEncoder'} + return data + + def __GPSL1TwoBitsEncoder_ToMap(self, obj): + data = {'type': 'GPSL1TwoBitsEncoder'} + return data + + def __GPSL2TwoBitsEncoder_ToMap(self, obj): + data = {'type': 'GPSL2TwoBitsEncoder'} + return data + + def __GPSL1L2TwoBitsEncoder_ToMap(self, obj): + data = {'type': 'GPSL1L2TwoBitsEncoder'} + return data + + def __MapTo_GPSL1BitEncoder(self, data): + return GPSL1BitEncoder() + + def __MapTo_GPSL2BitEncoder(self, data): + return GPSL2BitEncoder() + + def __MapTo_GPSL1L2BitEncoder(self, data): + return GPSL1L2BitEncoder() + + def __MapTo_GPSL1TwoBitsEncoder(self, data): + return GPSL1TwoBitsEncoder() + + def __MapTo_GPSL2TwoBitsEncoder(self, data): + return GPSL2TwoBitsEncoder() + + def __MapTo_GPSL1L2TwoBitsEncoder(self, data): + return GPSL1L2TwoBitsEncoder() + +factoryObject = ObjectFactory() diff --git a/peregrine/iqgen/bits/encoder_gps.py b/peregrine/iqgen/bits/encoder_gps.py new file mode 100644 index 0000000..99b46e1 --- /dev/null +++ b/peregrine/iqgen/bits/encoder_gps.py @@ -0,0 +1,194 @@ +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Valeri Atamaniouk +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + +""" +The :mod:`peregrine.iqgen.bits.encoder_gps` module contains classes and +functions related to generating GPS signal output. + +""" + +from peregrine.iqgen.bits.encoder_base import Encoder +from peregrine.iqgen.bits.encoder_1bit import BandBitEncoder +from peregrine.iqgen.bits.encoder_2bits import BandTwoBitsEncoder + + +class GPSL1BitEncoder(BandBitEncoder): + ''' + Generic single bit encoder for GPS L1 C/A signal + ''' + + def __init__(self, outputConfig): + ''' + Constructs GPS L1 C/A band single bit encoder object. + + Parameters + ---------- + outputConfig : object + Output parameters object. + ''' + super(GPSL1BitEncoder, self).__init__(outputConfig.GPS.L1.INDEX) + + +class GPSL2BitEncoder(BandBitEncoder): + ''' + Generic single bit encoder for GPS L2 Civil signal + ''' + + def __init__(self, outputConfig): + ''' + Constructs GPS L2 C band single bit encoder object. + + Parameters + ---------- + outputConfig : object + Output parameters object. + ''' + super(GPSL2BitEncoder, self).__init__(outputConfig.GPS.L2.INDEX) + + +class GPSL1L2BitEncoder(Encoder): + ''' + Generic single bit encoder for GPS L1 C/A and L2 Civil signals + ''' + + def __init__(self, outputConfig): + ''' + Constructs GPS L1 C/A and L2 C dual band single bit encoder object. + + Parameters + ---------- + outputConfig : object + Output parameters object. + ''' + super(GPSL1L2BitEncoder, self).__init__() + self.l1Index = outputConfig.GPS.L1.INDEX + self.l2Index = outputConfig.GPS.L2.INDEX + + def addSamples(self, sample_array): + ''' + Extracts samples of the supported band and coverts them into bit stream. + + Parameters + ---------- + sample_array : numpy.ndarray((4, N)) + Sample vectors ordered by band index. + + Returns + ------- + numpy.ndarray(dtype=numpy.uint8) + Array of type uint8 containing the encoded data. + ''' + band1_bits = BandBitEncoder.convertBand(sample_array[self.l1Index]) + band2_bits = BandBitEncoder.convertBand(sample_array[self.l2Index]) + n_samples = len(band1_bits) + + self.ensureExtraCapacity(n_samples * 2) + start = self.n_bits + end = start + 2 * n_samples + + self.bits[start + 0:end:2] = band1_bits + self.bits[start + 1:end:2] = band2_bits + self.n_bits = end + + if (self.n_bits >= Encoder.BLOCK_SIZE): + return self.encodeValues() + else: + return Encoder.EMPTY_RESULT + + +class GPSL1TwoBitsEncoder(BandTwoBitsEncoder): + ''' + Generic single bit encoder for GPS L1 C/A signal + ''' + + def __init__(self, outputConfig): + ''' + Constructs GPS L1 C/A band single bit encoder object. + + Parameters + ---------- + outputConfig : object + Output parameters object. + ''' + super(GPSL1TwoBitsEncoder, self).__init__(outputConfig.GPS.L1.INDEX) + + +class GPSL2TwoBitsEncoder(BandTwoBitsEncoder): + ''' + Generic single bit encoder for GPS L2 Civil signal + ''' + + def __init__(self, outputConfig): + ''' + Constructs GPS L2 C band single bit encoder object. + + Parameters + ---------- + outputConfig : object + Output parameters object. + ''' + super(GPSL2TwoBitsEncoder, self).__init__(outputConfig.GPS.L2.INDEX) + + +class GPSL1L2TwoBitsEncoder(Encoder): + ''' + Generic single bit encoder for GPS L1 C/A and L2 Civil signals + ''' + + def __init__(self, outputConfig): + ''' + Constructs GPS L1 C/A and L2 C dual band single bit encoder object. + + Parameters + ---------- + outputConfig : object + Output parameters object. + ''' + super(GPSL1L2TwoBitsEncoder, self).__init__() + self.l1Index = outputConfig.GPS.L1.INDEX + self.l2Index = outputConfig.GPS.L2.INDEX + + def addSamples(self, sample_array): + ''' + Extracts samples of the supported band and coverts them into bit stream. + + Parameters + ---------- + sample_array : numpy.ndarray((4, N)) + Sample vectors ordered by band index. + + Returns + ------- + numpy.ndarray(dtype=numpy.uint8) + Array of type uint8 containing the encoded data. + ''' + band1_samples = sample_array[self.l1Index] + band2_samples = sample_array[self.l2Index] + n_samples = len(band1_samples) + + # Signal signs and amplitude + signs1, amps1 = BandTwoBitsEncoder.convertBand(band1_samples) + signs2, amps2 = BandTwoBitsEncoder.convertBand(band2_samples) + + self.ensureExtraCapacity(n_samples * 4) + + bits = self.bits + start = self.n_bits + end = start + 4 * n_samples + bits[start + 0:end:4] = signs1 + bits[start + 1:end:4] = amps1 + bits[start + 2:end:4] = signs2 + bits[start + 3:end:4] = amps2 + self.n_bits = end + + if (self.n_bits >= Encoder.BLOCK_SIZE): + return self.encodeValues() + else: + return Encoder.EMPTY_RESULT diff --git a/peregrine/iqgen/bits/filter_bandpass.py b/peregrine/iqgen/bits/filter_bandpass.py new file mode 100644 index 0000000..753f6b8 --- /dev/null +++ b/peregrine/iqgen/bits/filter_bandpass.py @@ -0,0 +1,70 @@ +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Valeri Atamaniouk +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + +""" +The :mod:`peregrine.iqgen.bits.low_pass_filter` module contains classes and +functions related to generated signal attenuation. + +""" + +from scipy.signal.signaltools import lfiltic +from scipy.signal import cheby2, cheb2ord +from peregrine.iqgen.bits.filter_base import FilterBase + + +class BandPassFilter(FilterBase): + ''' + Chebyshev type 2 band-pass filter. + ''' + + def __init__(self, outputConfig, frequency_hz, bw_hz=1e6): + ''' + Initialize filter object. + + Parameters + ---------- + outputConfig : object + Output configuration parameters object + frequency_hz : float + Intermediate frequency in hertz + bw_hz : float, optional + Noise bandwidth in hertz + ''' + super(BandPassFilter, self).__init__(3., 40.) + + self.bw_hz = bw_hz + self.frequency_hz = frequency_hz + passBand_hz = bw_hz * 0.5 / outputConfig.SAMPLE_RATE_HZ + stopBand_hz = bw_hz * 0.6 / outputConfig.SAMPLE_RATE_HZ + mult = 2. / outputConfig.SAMPLE_RATE_HZ + order, wn = cheb2ord(wp=[(frequency_hz - passBand_hz) * mult, + (frequency_hz + passBand_hz) * mult], + ws=[(frequency_hz - stopBand_hz) * mult, + (frequency_hz + stopBand_hz) * mult], + gpass=self.passBandAtt_dbhz, + gstop=self.stopBandAtt_dbhz, + analog=False) + + b, a = cheby2(order, # Order of the filter + # Minimum attenuation required in the stop band in dB + self.stopBandAtt_dbhz, + wn, + btype="bandpass", + analog=False, + output='ba') + + self.a = a + self.b = b + self.zi = lfiltic(self.b, self.a, []) + + def __str__(self, *args, **kwargs): + return "BandPassFilter(center=%f, bw=%f, pb=%f, sp=%f)" % \ + (self.frequency_hz, self.bw_hz, + self.passBandAtt_dbhz, self.stopBandAtt_dbhz) diff --git a/peregrine/iqgen/bits/filter_base.py b/peregrine/iqgen/bits/filter_base.py new file mode 100644 index 0000000..179ddd4 --- /dev/null +++ b/peregrine/iqgen/bits/filter_base.py @@ -0,0 +1,71 @@ +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Valeri Atamaniouk +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + +""" +The :mod:`peregrine.iqgen.bits.filter_base` module contains classes and +functions related to generated signal attenuation. + +""" + +from scipy.signal import lfilter + + +class FilterBase(object): + + def __init__(self, passBandAtt_dbhz, stopBandAtt_dbhz): + ''' + Parameters + ---------- + passBandAtt_dbhz : float + Pass band attenutation in dB*Hz + stopBandAtt_dbhz: float + Stop band attenutation in dB*Hz + ''' + self.passBandAtt_dbhz = passBandAtt_dbhz + self.stopBandAtt_dbhz = stopBandAtt_dbhz + self.a = None + self.b = None + self.zi = None + + def getPassBandAtt(self): + ''' + Returns + ------- + float + Pass band attenuation in dB*Hz. + ''' + return self.passBandAtt_dbhz + + def getStopBandAtt(self): + ''' + Returns + ------- + float + Pass band attenuation in dB*Hz. + ''' + return self.stopBandAtt_dbhz + + def filter(self, data): + ''' + Performs noise reduction using Chebyshev type 2 IIR filter. + + Parameters + ---------- + data : array-like + Data samples before LPF processing + + Returns + ------- + array-like + Data samples after LPF processing + ''' + data_out, zo = lfilter(self.b, self.a, data, zi=self.zi) + self.zi = zo + return data_out diff --git a/peregrine/iqgen/bits/filter_lowpass.py b/peregrine/iqgen/bits/filter_lowpass.py new file mode 100644 index 0000000..37c6f5d --- /dev/null +++ b/peregrine/iqgen/bits/filter_lowpass.py @@ -0,0 +1,81 @@ +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Valeri Atamaniouk +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + +""" +The :mod:`peregrine.iqgen.bits.filter_lowpass` module contains classes and +functions related to generated signal attenuation. + +""" + +from scipy.signal.signaltools import lfiltic +from scipy.signal import cheby2, cheb2ord + +import logging + +from peregrine.iqgen.bits.filter_base import FilterBase + +logger = logging.getLogger(__name__) + + +class LowPassFilter(FilterBase): + ''' + Chebyshev type 2 low-pass filter. + + The filter simulates receiver lowpass effect: + - ~2MHz lowpass at 5*1.023MHz sampling pkg load signal + - b, a = cheby2(5, 40, 3e6/sample_freq) + + For 5.115MHz the coefficents are as follows: + - b = [0.082680, 0.245072, 0.397168, 0.397168 , 0.245072, 0.082680] + - a = [1.0000000, -0.3474596, 0.7770501, -0.0737540, 0.0922819, 0.0017200] + + ''' + + def __init__(self, outputConfig, frequency_hz=0.): + ''' + Initialize filter object. + + Parameters + ---------- + outputConfig : object + Output configuration parameters object + frequency_hz : float + Intermediate frequency + ''' + super(LowPassFilter, self).__init__(3., 40.) + + self.bw_hz = 1e3 + passBand_hz = self.bw_hz / outputConfig.SAMPLE_RATE_HZ + stopBand_hz = self.bw_hz * 1.1 / outputConfig.SAMPLE_RATE_HZ + mult = 2. / outputConfig.SAMPLE_RATE_HZ + order, wn = cheb2ord(wp=passBand_hz * mult, + ws=stopBand_hz * mult, + gpass=self.passBandAtt_dbhz, + gstop=self.stopBandAtt_dbhz, + analog=False) + self.order = order + self.wn = wn + + b, a = cheby2(order, # Order of the filter + # Minimum attenuation required in the stop band in dB + self.stopBandAtt_dbhz, + wn, + btype="lowpass", + analog=False, + output='ba') + + self.a = a + self.b = b + self.zi = lfiltic(self.b, self.a, []) + + def __str__(self, *args, **kwargs): + return "LowPassFilter(bw=%f, pb=%f, sp=%f, order=%d, wn=%s)" % \ + (self.bw_hz, self.passBandAtt_dbhz, + self.stopBandAtt_dbhz, self.order, str(self.wn)) diff --git a/peregrine/iqgen/bits/message_block.py b/peregrine/iqgen/bits/message_block.py new file mode 100644 index 0000000..425efff --- /dev/null +++ b/peregrine/iqgen/bits/message_block.py @@ -0,0 +1,85 @@ +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Valeri Atamaniouk +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + +""" +The :mod:`peregrine.iqgen.bits.message_block` module contains classes and +functions related to providing predefined symbol contents. + +""" +import numpy + + +class Message(object): + ''' + Message that is a block of bits + ''' + + def __init__(self, messageData): + ''' + Constructs message object. + + Parameters + ---------- + messageData : array-like + Array with message bits. Bit 0 is encoded with 1, bit 1 is encoded with -1 + ''' + super(Message, self).__init__() + self.messageData = messageData[:] + self.messageLen = len(self.messageData) + tmp = numpy.asarray(self.messageData, dtype=numpy.uint8) + tmp *= -2 + tmp -= 1 + self.bits = tmp.astype(numpy.uint8) + + def __str__(self, *args, **kwargs): + ''' + Formats object as string literal + + Returns + ------- + string + String representation of the object + ''' + return "Block: length=%d" % len(self.bits) + + def getDataBits(self, dataAll_idx): + ''' + Generates vector of data bits corresponding to input index + + Parameters + ---------- + dataAll_idx : numpy.ndarray(dtype=numpy.int64) + Vector of bit indexes + + Returns + ------- + numpy.ndarray(dtype=numpy.uint8) + Vector of data bits + ''' + # numpy.take degrades performance a lot over time. + # return numpy.take(self.bits, dataAll_idx , mode='wrap') + return self.bits[dataAll_idx % self.messageLen] + + def getBit(self, bitIndex): + ''' + Provides bit at a given index + + Parameters + ---------- + bitIndex : long + Bit index + + Returns + ------- + int + Bit value: 1 for bit 0 and -1 for bit 1 + ''' + + return self.messageData[bitIndex % self.messageLen] diff --git a/peregrine/iqgen/bits/message_cnav.py b/peregrine/iqgen/bits/message_cnav.py new file mode 100644 index 0000000..c280500 --- /dev/null +++ b/peregrine/iqgen/bits/message_cnav.py @@ -0,0 +1,222 @@ +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Valeri Atamaniouk +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + +''' +The :mod:`peregrine.iqgen.bits.message_cnav` module contains classes and +functions related to generating stub GPS CNAV messages. +''' + +import numpy +from swiftnav.cnav_msg import CNavRawMsg +import logging +G1 = 0171 # generator polinomial for p1 +G2 = 0133 # generator polinomial for p2 + +logger = logging.getLogger(__name__) + + +def generate27Vector(g1, g2): + ''' + Helper method for convolution encoder lookup table generation. + + Parameters + ---------- + g1 : int + First polynomial coefficient + g2 : int + Second polynomial coefficient + + Results + ------- + numpy.ndvector(shape=(128,2),dtype=numpy.uint8) + Lookup matrix for convolution encoder + ''' + + def parity6(value): + ''' + Helper for computing parity of 6-bit value. + + Parameters + ---------- + value : int + 6-bit integer value + + Results + ------- + int + Parity bit: 0 or 1. + ''' + return (0x6996 >> ((value ^ (value >> 4)) & 15)) & 1 + + vectorG = numpy.ndarray((128, 2), dtype=numpy.uint8) + for i in range(128): + vectorG[i][0] = parity6(i & g1) + vectorG[i][1] = parity6(i & g2) + + return vectorG + + +class ConvEncoder27(object): + ''' + Convolution encoder class. + + Standard 2-7 convolution encoder implementation. + ''' + + DEFAULT_VECTOR_G = generate27Vector(G1, G2) + + def __init__(self, g1=G1, g2=G2, state=0): + self.g1 = g1 + self.g2 = g2 + self.state = state + vectorG = ConvEncoder27.DEFAULT_VECTOR_G if g1 == G1 and g2 == G2 \ + else generate27Vector(g1, g2) + self.vectorG = vectorG + + def encode(self, bitArray): + ''' + Encodes source bit array. + + This method updates the encoder state during processing. + + Parameters + ---------- + bitArray : array-like + Array of bit values. Can be integers or booleans. + Returns + ------- + numpy.ndarray(shape(len(bitArray)), dtype=numpy.uint8) + Encoded output + ''' + result = numpy.ndarray((len(bitArray) * 2), dtype=numpy.uint8) + state = self.state + dstIndex = 0 + vectorG = self.vectorG + + for srcBit in bitArray: + state = (srcBit << 6) | (state >> 1) + result[dstIndex:dstIndex + 2] = vectorG[state] + dstIndex += 2 + + self.state = state + return result + + +class Message(object): + ''' + GPS LNAV message block. + + The object provides proper-formatted CNAV messages with random contents. + ''' + + def __init__(self, prn, tow0=2, n_msg=0, n_prefixBits=50): + ''' + Constructs message object. + + Parameters + ---------- + prn : int + Satellite PRN + tow0 : int + Time of week in 6-second units for the first message + n_msg : int, optional + Number of messages to generate for output + n_prefixBits : int, optional + Number of bits to issue before the first message + ''' + super(Message, self).__init__() + + if tow0 & 1: + logger.error("Initial ToW is not multiple of 2") + + self.prn = prn + self.tow0 = tow0 + self.n_msg0 = n_msg + self.n_prefixBits = n_prefixBits + + self.encoder = ConvEncoder27() + self.msgCount = 0 + self.messageLen = n_prefixBits * 2 + self.symbolData = numpy.zeros(self.messageLen, dtype=numpy.uint8) + + prefixBits = numpy.zeros(self.n_prefixBits, dtype=numpy.uint8) + prefixBits[0::2] = 1 + self.symbolData[:] = self.encoder.encode(prefixBits) + self.nextTow = tow0 + self.addMessages(n_msg) + + def __str__(self, *args, **kwargs): + ''' + Formats object as string literal + + Returns + ------- + string + String representation of the object + ''' + return "GPS CNAV: prn=%d pref=%d tow=%d" % \ + (self.prn, self.n_prefixBits, self.nextTow) + + def getDataBits(self, dataAll_idx): + ''' + Generates vector of data bits corresponding to input index + + Parameters + ---------- + dataAll_idx : numpy.ndarray(dtype=numpy.int64) + Vector of bit indexes + + Returns + ------- + numpy.ndarray(dtype=numpy.uint8) + Vector of data bits + ''' + + lastIdx = dataAll_idx[-1] + if lastIdx >= self.messageLen: + # Grow data bits + delta = lastIdx - self.messageLen + 1 + newMsgCount = delta / 600 + if delta % 600: + newMsgCount += 1 + self.addMessages(newMsgCount) + + # numpy.take degrades performance a lot over time. + # return numpy.take(self.symbolData, dataAll_idx , mode='wrap') + return self.symbolData[dataAll_idx] + + def addMessages(self, newMsgCount): + ''' + Generate additional CNAV messages + + This method generates and encodes additional CNAV messages. The message + contents is encoded using 2-7 convolution encoder and added to the internal + buffer. + + Parameters + ---------- + newMsgCount : int + Number of messages to generate + ''' + newMessageLen = newMsgCount * 600 + self.messageLen + newSymbolData = numpy.ndarray(newMessageLen, dtype=numpy.uint8) + newSymbolData[:self.messageLen] = self.symbolData + for i in range(self.messageLen, newMessageLen, 600): + logger.info("Generating CNAV message: prn=%d tow=%d msg_id=%d" % + (self.prn, self.nextTow, 0)) + cnav_msg = CNavRawMsg.generate(self.prn, 0, self.nextTow) + self.nextTow += 2 + if self.nextTow == 7 * 24 * 60 * 10: + self.nextTow = 0 + encoded = self.encoder.encode(cnav_msg) + newSymbolData[i:i + 600] = encoded + self.messageLen = newMessageLen + self.symbolData = newSymbolData + self.msgCount += newMsgCount diff --git a/peregrine/iqgen/bits/message_const.py b/peregrine/iqgen/bits/message_const.py new file mode 100644 index 0000000..dc3e135 --- /dev/null +++ b/peregrine/iqgen/bits/message_const.py @@ -0,0 +1,65 @@ +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Valeri Atamaniouk +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + +import numpy + +""" +The :mod:`peregrine.iqgen.bits.message_const` module contains classes and +functions related to non-changing symbol contents. + +""" + + +class Message(object): + ''' + Message consisting of same bits + ''' + + def __init__(self, bitValue): + ''' + Initializes object. + + Parameters + ---------- + bitValue : int + Value for the bits. 1 for 0 bits, -1 for 1 bits. + ''' + super(Message, self).__init__() + self.bitValue = bitValue + self.binValue = 1 if bitValue < 0 else 0 + + def __str__(self, *args, **kwargs): + ''' + Formats object as string literal + + Returns + ------- + string + String representation of the object + ''' + return "Const: bit value=%d" % self.binValue + + def getDataBits(self, dataAll_idx): + ''' + Generates vector of data bits corresponding to input index + + Parameters + ---------- + dataAll_idx : numpy.ndarray(dtype=numpy.int64) + Vector of bit indexes + + Returns + ------- + numpy.ndarray(dtype=numpy.uint8) + Vector of data bits + ''' + result = numpy.ndarray(len(dataAll_idx), dtype=numpy.uint8) + result.fill(self.binValue) + return result diff --git a/peregrine/iqgen/bits/message_factory.py b/peregrine/iqgen/bits/message_factory.py new file mode 100644 index 0000000..67c5034 --- /dev/null +++ b/peregrine/iqgen/bits/message_factory.py @@ -0,0 +1,124 @@ +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Valeri Atamaniouk +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + + +""" +The :mod:`peregrine.iqgen.bits.messager_factory` module contains classes and +functions related to object factory for message objects. + +""" + +from peregrine.iqgen.bits.message_block import Message as BlockMessage +from peregrine.iqgen.bits.message_cnav import Message as CNAVMessage +from peregrine.iqgen.bits.message_const import Message as ConstMessage +from peregrine.iqgen.bits.message_lnav import Message as LNAVMessage +from peregrine.iqgen.bits.message_zeroone import Message as ZeroOneMessage + + +class ObjectFactory(object): + ''' + Object factory for message objects. + ''' + + def __init__(self): + super(ObjectFactory, self).__init__() + + def toMapForm(self, obj): + t = type(obj) + if t is BlockMessage: + return self.__BlockMessage_ToMap(obj) + elif t is CNAVMessage: + return self.__CNAVMessage_ToMap(obj) + elif t is ConstMessage: + return self.__ConstMessage_ToMap(obj) + elif t is LNAVMessage: + return self.__LNAVMessage_ToMap(obj) + elif t is ZeroOneMessage: + return self.__ZeroOneMessage_ToMap(obj) + else: + raise ValueError("Invalid object type") + + def fromMapForm(self, data): + t = data['type'] + if t == 'BlockMessage': + return self.__MapTo_BlockMessage(data) + elif t == 'CNAVMessage': + return self.__MapTo_CNAVMessage(data) + elif t == 'ConstMessage': + return self.__MapTo_ConstMessage(data) + elif t == 'LNAVMessage': + return self.__MapTo_LNAVMessage(data) + elif t == 'ZeroOneMessage': + return self.__MapTo_ZeroOneMessage(data) + else: + raise ValueError("Invalid object type") + + def __BlockMessage_ToMap(self, obj): + data = {'type': 'BlockMessage', + 'data': obj.messageData} + return data + + def __CNAVMessage_ToMap(self, obj): + data = {'type': 'CNAVMessage', + 'prn': obj.prn, + 'n_prefixBits': obj.n_prefixBits, + 'n_msg0': obj.n_msg0, + 'tow0': obj.tow0} + return data + + def __ConstMessage_ToMap(self, obj): + data = {'type': 'ConstMessage', + 'bitValue': obj.bitValue} + return data + + def __LNAVMessage_ToMap(self, obj): + data = {'type': 'LNAVMessage', + 'prn': obj.prn, + 'n_prefixBits': obj.n_prefixBits, + 'n_msg0': obj.n_msg0, + 'tow0': obj.tow0} + return data + + def __ZeroOneMessage_ToMap(self, obj): + data = {'type': 'ZeroOneMessage'} + return data + + def __MapTo_BlockMessage(self, data): + messageData = data['data'] + return BlockMessage(messageData) + + def __MapTo_CNAVMessage(self, data): + prn = data['prn'] + n_prefixBits = data['n_prefixBits'] + n_msg0 = data['n_msg0'] + tow0 = data['tow0'] + return CNAVMessage(prn=prn, + tow0=tow0, + n_msg=n_msg0, + n_prefixBits=n_prefixBits) + + def __MapTo_ConstMessage(self, data): + bitValue = data['bitValue'] + return ConstMessage(bitValue) + + def __MapTo_LNAVMessage(self, data): + prn = data['prn'] + n_prefixBits = data['n_prefixBits'] + n_msg0 = data['n_msg0'] + tow0 = data['tow0'] + return LNAVMessage(prn=prn, + tow0=tow0, + n_msg=n_msg0, + n_prefixBits=n_prefixBits) + + def __MapTo_ZeroOneMessage(self, data): + return ZeroOneMessage() + +factoryObject = ObjectFactory() diff --git a/peregrine/iqgen/bits/message_lnav.py b/peregrine/iqgen/bits/message_lnav.py new file mode 100644 index 0000000..6a595ed --- /dev/null +++ b/peregrine/iqgen/bits/message_lnav.py @@ -0,0 +1,273 @@ +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Valeri Atamaniouk +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + +''' +The :mod:`peregrine.iqgen.bits.message_lnav` module contains classes and +functions related to generating stub GPS LNAV messages. +''' + +import numpy +import logging +from swiftnav.bits import parity + +logger = logging.getLogger(__name__) + + +class Message(object): + ''' + GPS LNAV message generator + ''' + + def __init__(self, prn, tow0=1, n_msg=0, n_prefixBits=50): + ''' + Constructs message object. + + Parameters + ---------- + prn : int + Satellite PRN + tow0 : int + Time of week in 6-second units for the first message + n_msg : int, optional + Number of messages to generate for output + n_prefixBits : int, optional + Number of bits to issue before the first message + ''' + super(Message, self).__init__() + self.prn = prn + self.n_prefixBits = n_prefixBits + self.n_msg0 = n_msg + self.tow0 = tow0 + self.messageCount = 0 + self.messageLen = n_prefixBits + self.nextTow = tow0 + self.nextMsgId = 0 + self.messageBits = numpy.zeros(n_prefixBits, dtype=numpy.uint8) + words = (n_prefixBits + 29) / 30 + if words: + tmp = numpy.zeros(words * 30, dtype=numpy.uint8) + tmp[1::2] = 1 + if words > 1: + self.updateParity(tmp[0:30]) + for i in range(1, words - 1): + self.updateParity(tmp[i * 30 - 2: i * 30 + 30]) + self.updateParity(tmp[words * 30 - 32: words * 30], True) + else: + self.updateParity(tmp[0: 30], True) + self.messageBits[:] = tmp[-n_prefixBits:] + self.msgCount = 0 + self.a8 = numpy.ndarray(1, dtype=numpy.uint8) + self.a32 = numpy.ndarray(1, dtype=numpy.dtype('>u4')) + + def __str__(self, *args, **kwargs): + ''' + Formats object as string literal + + Returns + ------- + string + String representation of the object + ''' + return "GPS LNAV: prn=%d pref=%d tow=%d" % \ + (self.prn, self.n_prefixBits, self.nextTow) + + def getDataBits(self, dataAll_idx): + ''' + Generates vector of data bits corresponding to input index + + Parameters + ---------- + dataAll_idx : numpy.ndarray(dtype=numpy.int64) + Vector of bit indexes + + Returns + ------- + numpy.ndarray(dtype=numpy.uint8) + Vector of data bits + ''' + + lastIdx = dataAll_idx[-1] + if lastIdx >= self.messageLen: + # Grow data bits + delta = lastIdx - self.messageLen + 1 + newMsgCount = delta / 300 + if delta % 300: + newMsgCount += 1 + self.addMessages(newMsgCount) + + # numpy.take degrades performance a lot over time. + # return numpy.take(self.symbolData, dataAll_idx , mode='wrap') + return self.messageBits[dataAll_idx] + + def addMessages(self, newMsgCount): + ''' + Generate additional LNAV messages + + This method generates and encodes additional LNAV messages. The message + contents is added to the internal buffer. + + Parameters + ---------- + newMsgCount : int + Number of messages to generate + ''' + newMessageLen = newMsgCount * 300 + self.messageLen + newMessageData = numpy.ndarray(newMessageLen, dtype=numpy.uint8) + newMessageData[:self.messageLen] = self.messageBits + for i in range(self.messageLen, newMessageLen, 300): + logger.info("Generating LNAV message: prn=%d tow=%d msg_id=%d" % + (self.prn, self.nextTow, self.nextMsgId)) + lnav_msg = self.generateLNavMessage() + newMessageData[i:i + 300] = lnav_msg + self.messageLen = newMessageLen + self.messageBits = newMessageData + self.msgCount += newMsgCount + + def generateLNavMessage(self): + ''' + Produces additional GPS LNAV message. + + Returns + ------- + numpy.ndarray(shape=300, dtype=numpy.uint8) + Message bits. + ''' + msgData = numpy.zeros(300, dtype=numpy.uint8) + msgData[1::2] = 1 # Zero + one everywhere + + # TLM word + self.fillTlmWord(msgData[0:30], 0) + self.updateParity(msgData[0:30]) + # logger.debug("TLM: %s" % msgData[0:30]) + + # TOW word + self.fillTowWord(msgData[30:60], self.nextTow) + self.nextTow += 1 + if self.nextTow == 7 * 24 * 60 * 10: + self.nextTow = 0 + self.updateParity(msgData[28:60], True) + # logger.debug("TOW: %s" % msgData[30:60]) + + self.updateParity(msgData[58:90]) + self.updateParity(msgData[88:120]) + self.updateParity(msgData[118:150]) + self.updateParity(msgData[148:180]) + self.updateParity(msgData[178:210]) + self.updateParity(msgData[208:240]) + self.updateParity(msgData[238:270]) + self.updateParity(msgData[268:300], True) + + return msgData + + def getBits(self, value, nBits): + ''' + Converts integer into bit array + + Parameters + ---------- + value : int + Integer value + nBits : number of bits to produce + + Returns + ------- + numpy.ndarray(shape=(`nBits`), dtype=numpy.uint8) + Parameter `value` represented as a bit array + ''' + if nBits <= 8: + self.a8[0] = value + result = numpy.unpackbits(self.a8) + else: + self.a32[0] = value + result = numpy.unpackbits(self.a32.view(dtype=numpy.uint8)) + return result[-nBits:] + + def fillTlmWord(self, wordBits, msgId=0): + ''' + Fills in TLM word contents. + + Parameters + ---------- + wordBits : numpy.ndarray(shape=30, type=numpy.uint8) + Destination array + ''' + wordBits[0:8] = self.getBits(0b10001011, 8) # Preamble + wordBits[8:22] = self.getBits(msgId, 14) # TML message + wordBits[22] = 0 # Reserved + wordBits[23] = 0 # Integrity + return + + def fillTowWord(self, wordBits, tow): + ''' + Fills in TOW word contents. + + Parameters + ---------- + wordBits : numpy.ndarray(shape=30, type=numpy.uint8) + Destination array + ''' + wordBits[0:17] = self.getBits(tow, 17) # TOW count in 6 second units + wordBits[17] = 0 # Alert Flag + wordBits[18] = 0 # Anti-Spoof flag + wordBits[19:22] = self.getBits(0, 3) # Sub-frame ID + return + + def updateParity(self, dataBits, resolve=False): + ''' + Updates data bits and computes parity. + + When 32 bits are provided, they are used for parity computation and for + bit inversion. + + Parameters + ---------- + dataBits : numpy.ndarray(dtype=numpy.uint8) + 30 or 32 element array + resolve: bool, optional + When specified, bits d23 and d24 of the GPS word are updated to ensure + that parity bits d29 and d30 are zeros. + ''' + packed = numpy.packbits(dataBits) + acc = (packed[0] << 24) | (packed[1] << 16) | \ + (packed[2] << 8) | packed[3] + if len(dataBits) == 30: + acc >>= 2 + elif acc & 0x40000000: + acc ^= 0x3FFFFFC0 + dataBits[-30:-6] ^= 1 + + # D29 = D30*^d1^d3^d5^d6^d7^d9^d10^d14^d15^d16^d17^d18^d21^d22^d24 + d29 = parity(acc & 0b01101011101100011111001101000000) + # D30 = D29*^d3^d5^d6^d8^d9^d10^d11^d13^d15^d19^d22^d23^d24 + d30 = parity(acc & 0b10001011011110101000100111000000) + + if resolve: + if d29: + acc ^= 0x80 + d29 = False + d30 = not d30 + dataBits[-8] ^= 1 + if d30: + acc ^= 0x40 + d30 = False + dataBits[-7] ^= 1 + + # D25 = D29*^d1^d2^d3^d5^d6^d10^d11^d12^d13^d14^d17^d18^d20^d23 + dataBits[-6] = parity(acc & 0b10111011000111110011010010000000) + # D26 = D30*^d2^d3^d4^d6^d7^d11^d12^d13^d14^d15^d18^d19^d21^d24 + dataBits[-5] = parity(acc & 0b01011101100011111001101001000000) + # D27 = D29*^d1^d3^d4^d5^d7^d8^d12^d13^d14^d15^d16^d19^d20^d22 + dataBits[-4] = parity(acc & 0b10101110110001111100111000000000) + # D28 = D30*^d2^d4^d5^d6^d8^d9^d13^d14^d15^d16^d17^d20^d21^d23 + dataBits[-3] = parity(acc & 0b01010111011000111110011010000000) + # D29 = D30*^d1^d3^d5^d6^d7^d9^d10^d14^d15^d16^d17^d18^d21^d22^d24 + dataBits[-2] = d29 + # D30 = D29*^d3^d5^d6^d8^d9^d10^d11^d13^d15^d19^d22^d23^d24 + dataBits[-1] = d30 diff --git a/peregrine/iqgen/bits/message_zeroone.py b/peregrine/iqgen/bits/message_zeroone.py new file mode 100644 index 0000000..0882f53 --- /dev/null +++ b/peregrine/iqgen/bits/message_zeroone.py @@ -0,0 +1,59 @@ +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Valeri Atamaniouk +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + +""" +The :mod:`peregrine.iqgen.bits.message_zeroone` module contains classes and +functions related to symbol contents that flips the value every other bit. + +""" + +import numpy + + +class Message(object): + ''' + Message that contains zeros and ones + ''' + + def __init__(self): + ''' + Constructs object. + ''' + super(Message, self).__init__() + self.bits = numpy.asarray([0, 1], dtype=numpy.uint8) + + def __str__(self, *args, **kwargs): + ''' + Formats object as string literal + + Returns + ------- + string + String representation of the object + ''' + return "ZeroOne" + + def getDataBits(self, dataAll_idx): + ''' + Generates vector of data bits corresponding to input index + + Parameters + ---------- + dataAll_idx : numpy.ndarray(dtype=numpy.int64) + Vector of bit indexes + + Returns + ------- + numpy.ndarray(dtype=numpy.uint8) + Vector of data bits + ''' + # numpy.take degrades performance a lot over time. + # return numpy.take(self.bits, dataAll_idx , mode='wrap') + return self.bits[dataAll_idx & 1] diff --git a/peregrine/iqgen/bits/prn_gps_l1ca.py b/peregrine/iqgen/bits/prn_gps_l1ca.py new file mode 100644 index 0000000..96893d6 --- /dev/null +++ b/peregrine/iqgen/bits/prn_gps_l1ca.py @@ -0,0 +1,86 @@ +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Valeri Atamaniouk +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + +import numpy + +""" +The :mod:`peregrine.iqgen.bits.prn_gps_l1ca` module contains classes and +functions related to GPS L1 C/A PRN processing + +""" + +import peregrine.include.generateCAcode + +caCodes = peregrine.include.generateCAcode.caCodes + + +class PrnCode(object): + ''' + GPS L1 C/A code object + ''' + CODE_LENGTH = 1023 + CODE_FREQUENCY_HZ = 1023e3 + + def __init__(self, prnNo): + ''' + Initializes object. + + Parameters + ---------- + prnNo : int + SV identifier + ''' + super(PrnCode, self).__init__() + self.caCode = caCodes[prnNo - 1][:] + tmp = numpy.asarray(self.caCode, dtype=numpy.int8) + tmp -= 1 + tmp /= -2 + self.binCode = tmp + self.prnNo = prnNo + self.bitLookup = numpy.asarray([1, -1], dtype=numpy.int8) + + def getCodeBits(self, chipIndex_all): + ''' + Parameters + ---------- + chipIndex_all : numpy.ndarray(dtype=numpy.long) + Vector of chip indexes + + Returns + ------- + numpy.ndarray(dtype=numpy.uint8) + Vector of code chip bits + ''' + # numpy.take degrades performance a lot over time. + # return numpy.take(self.binCode, chipIndex_all, mode='wrap') + return self.binCode[chipIndex_all % len(self.binCode)] + + def combineData(self, chipIndex_all, dataBits): + ''' + Mixes in code chip and data + + Parameters + ---------- + chipIndex_all : numpy.ndarray(dtype=numpy.long) + Chip indexes + dataBits : numpy.ndarray(dtype=numpy.uint8) + Data bits + + Returns + ------- + numpy.ndarray(dtype=numpy.int8) + Vector of data bits modulated by chips + ''' + chipBits = self.getCodeBits(chipIndex_all) + combined = numpy.bitwise_xor(chipBits, dataBits) + # numpy.take degrades performance a lot over time. + # result = numpy.take(self.bitLookup, combined) + result = self.bitLookup[combined] + return result diff --git a/peregrine/iqgen/bits/prn_gps_l2c.py b/peregrine/iqgen/bits/prn_gps_l2c.py new file mode 100644 index 0000000..e424e58 --- /dev/null +++ b/peregrine/iqgen/bits/prn_gps_l2c.py @@ -0,0 +1,206 @@ +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Valeri Atamaniouk +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + + +""" +The :mod:`peregrine.iqgen.bits.prn_gps_l2c` module contains classes and +functions related to GPS L2C PRN processing + +""" + +import numpy + +from peregrine.include.generateL2CMcode import L2CMCodes + + +class PrnCode(object): + ''' + Combined GPS L2 CM and CL code object + ''' + + class CM_Code(object): + ''' + GPS L2 Civil Medium code object + ''' + CODE_LENGTH = 10230 + CODE_FREQUENCY_HZ = 511.5e3 + + def __init__(self, prnNo): + ''' + Initializes object. + + Parameters + ---------- + prnNo : int + SV identifier + ''' + super(PrnCode.CM_Code, self).__init__() + self.caCode = L2CMCodes[prnNo - 1][:] + self.binCode = numpy.ndarray( + PrnCode.CM_Code.CODE_LENGTH, dtype=numpy.bool) + self.binCode[:] = numpy.asarray(self.caCode) < 0 + self.prnNo = prnNo + + def getCodeBits(self): + return self.binCode + + def getCodeBit(self, codeBitIndex): + ''' + Returns chip value by index. + + Parameters + ---------- + chipIndex : long + Chip index + + Returns + ------- + int + Chip value by index + ''' + return self.caCode[codeBitIndex % self.CODE_LENGTH] + + class CL_Code(object): + ''' + GPS L2 Civil Long code object + ''' + CODE_LENGTH = 767250 + CODE_FREQUENCY_HZ = 511.5e3 + + def __init__(self, prnNo, codeType): + ''' + Initializes object. + + Parameters + ---------- + prnNo : int + SV identifier + codeType : string + Type of the code: '01', '1', '0' + ''' + super(PrnCode.CL_Code, self).__init__() + self.prnNo = prnNo + self.binCode = numpy.ndarray( + PrnCode.CL_Code.CODE_LENGTH, dtype=numpy.bool) + if codeType == '01': + self.binCode.fill(False) + self.binCode[1::2].fill(True) + elif codeType == '1': + self.binCode.fill(True) + elif codeType == '0': + self.binCode.fill(False) + else: + raise ValueError('Unsupported GPS L2 CL generator type %s ' % + str(codeType)) + + def getCodeBits(self): + return self.binCode + + def getCodeBit(self, codeBitIndex): + ''' + Returns chip value by index. + + Currently GPS L2 CL code can be pseudo-random + + Parameters + ---------- + chipIndex : long + Chip index + + Returns + ------- + int + Chip value by index + ''' + if (codeBitIndex & 1 != 0): + return -1 + else: + return 1 + + CODE_LENGTH = CL_Code.CODE_LENGTH * 2 + CODE_FREQUENCY_HZ = 1023e3 + + def __init__(self, prnNo, clCodeType): + ''' + Initializes object. + + Parameters + ---------- + prnNo : int + SV identifier + clCodeType : string + Type of the code: '01', '1', '0' + ''' + super(PrnCode, self).__init__() + self.cl = PrnCode.CL_Code(prnNo, clCodeType) + self.cm = PrnCode.CM_Code(prnNo) + self.bitLookup = numpy.asarray([1, -1], dtype=numpy.int8) + tmp = numpy.ndarray(PrnCode.CL_Code.CODE_LENGTH * 2, dtype=numpy.uint8) + tmp[1::2] = self.cl.getCodeBits() + for i in range(0, PrnCode.CL_Code.CODE_LENGTH * 2, PrnCode.CM_Code.CODE_LENGTH * 2): + tmp[i:i + PrnCode.CM_Code.CODE_LENGTH * 2:2] = self.cm.getCodeBits() + self.binCode = tmp + self.prnNo = prnNo + + def getCodeBits(self, chipIndex_all): + ''' + Parameters + ---------- + chipIndex_all : numpy.ndarray(dtype=numpy.long) + Vector of chip indexes + + Returns + ------- + numpy.ndarray(dtype=numpy.uint8) + Vector of code chip bits + ''' + # numpy.take degrades performance a lot over time. + # return numpy.take(self.binCode, chipIndex_all, mode='wrap') + return self.binCode[chipIndex_all % PrnCode.CODE_LENGTH] + + def combineData(self, chipIndex_all, dataBits): + ''' + Mixes in code chip and data + + Parameters + ---------- + chipIndex_all : numpy.ndarray(dtype=numpy.long) + Chip indexes + dataBits : numpy.ndarray(dtype=numpy.uint8) + Data bits + + Returns + ------- + numpy.ndarray(dtype=numpy.int8) + Vector of data bits modulated by chips + ''' + chipBits = self.getCodeBits(chipIndex_all) + tmp = dataBits.copy() + oddChips = chipIndex_all & 1 == 0 + # print "idx=", chipIndex_all + # print "odd=", oddChips + # print "TMP1", tmp + tmp = tmp & oddChips + # print "TMP2", tmp + combined = numpy.bitwise_xor(chipBits, tmp) + # numpy.take degrades performance a lot over time. + # result = numpy.take(self.bitLookup, combined) + result = self.bitLookup[combined] + return result + + def __getCodeBit(self, codeBitIndex): + ''' + For GPS L2C code bits are taken from CM and CL codes in turn. + ''' + idx = long(codeBitIndex) + if idx & 1 != 0: + return self.cl.getCodeBit(idx / 2) + else: + return self.cm.getCodeBit(idx / 2) diff --git a/peregrine/iqgen/bits/satellite_base.py b/peregrine/iqgen/bits/satellite_base.py new file mode 100644 index 0000000..e5bfa0f --- /dev/null +++ b/peregrine/iqgen/bits/satellite_base.py @@ -0,0 +1,138 @@ +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Valeri Atamaniouk +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + +""" +The :mod:`peregrine.iqgen.bits.satellite_base` module contains classes and +functions related to satellite base object. + +""" +from peregrine.iqgen.bits.doppler_poly import zeroDoppler +from peregrine.iqgen.bits.amplitude_poly import AmplitudePoly + + +class Satellite(object): + ''' + Satellite object. + + Satellite object combines speed/position computation and data generation for + all supported bands. + ''' + + def __init__(self, svName): + ''' + Constructor. + + Parameters + ---------- + svName : string + Satellite name + ''' + super(Satellite, self).__init__() + self.svName = svName + self.doppler = zeroDoppler(0., 0., 1.) + self.amplitude = AmplitudePoly(()) + + def getDoppler(self): + ''' + Returns doppler object. + + Returns + ------- + object + Doppler control object + ''' + return self.doppler + + def setDoppler(self, doppler): + ''' + Changes doppler object. + + Parameters + ------- + doppler : object + Doppler control object + ''' + self.doppler = doppler + + def getSvName(self): + ''' + Returns satellite name. + + Returns + ------- + string + Satellite name + ''' + return self.svName + + def setAmplitude(self, amplitude): + ''' + Changes amplitude + + Parameters + ---------- + amplitude : float + amplitude value for signal generation + ''' + self.amplitude = amplitude + + def getAmplitude(self): + ''' + Provides amplitude object + + Returns + ------- + object + Amplitude object + ''' + return self.amplitude + + def __str__(self): + return self.getSvName() + + def __repr__(self): + return self.getSvName() + + def getBatchSignals(self, userTimeAll_s, samples, outputConfig): + ''' + Generates signal samples. + + Parameters + ---------- + userTimeAll_s : numpy.ndarray(n_samples, dtype=numpy.float64) + Vector of observer's timestamps in seconds for the interval start. + samples : numpy.ndarray((4, n_samples)) + Array to which samples are added. + outputConfig : object + Output configuration object. + + Returns + ------- + list + Debug information + ''' + raise NotImplementedError() + + def isBandEnabled(self, bandIndex, outputConfig): + ''' + Checks if particular band is supported and enabled. + + Parameters + ---------- + bandIndex : int + Signal band index + outputConfig : object + Output configuration + + Returns: + bool + True, if the band is supported and enabled; False otherwise. + ''' + return False diff --git a/peregrine/iqgen/bits/satellite_factory.py b/peregrine/iqgen/bits/satellite_factory.py new file mode 100644 index 0000000..f62fc15 --- /dev/null +++ b/peregrine/iqgen/bits/satellite_factory.py @@ -0,0 +1,81 @@ +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Valeri Atamaniouk +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + + +""" +The :mod:`peregrine.iqgen.bits.satellite_factory` module contains classes and +functions related to object factory for satellite objects. + +""" + +from peregrine.iqgen.bits.satellite_gps import GPSSatellite +from peregrine.iqgen.bits.amplitude_factory import factoryObject as amplitudeOF +from peregrine.iqgen.bits.doppler_factory import factoryObject as dopplerOF +from peregrine.iqgen.bits.message_factory import factoryObject as messageOF + + +class ObjectFactory(object): + ''' + Object factory for satellite types. + ''' + + def __init__(self): + super(ObjectFactory, self).__init__() + + def toMapForm(self, obj): + t = type(obj) + if t is GPSSatellite: + return self.__GPSSatellite_ToMap(obj) + else: + raise ValueError("Invalid object type") + + def fromMapForm(self, data): + t = data['type'] + if t == 'GPSSatellite': + return self.__MapTo_GPSSatellite(data) + else: + raise ValueError("Invalid object type") + + def __GPSSatellite_ToMap(self, obj): + data = {'type': 'GPSSatellite', + 'prn': obj.prn, + 'amplitude': amplitudeOF.toMapForm(obj.getAmplitude()), + 'l1caEnabled': obj.isL1CAEnabled(), + 'l2cEnabled': obj.isL2CEnabled(), + 'l1caMessage': messageOF.toMapForm(obj.getL1CAMessage()), + 'l2cMessage': messageOF.toMapForm(obj.getL2CMessage()), + 'doppler': dopplerOF.toMapForm(obj.getDoppler()), + 'l2clCodeType': obj.getL2CLCodeType(), + 'codeDopplerIgnored': obj.isCodeDopplerIgnored() + } + return data + + def __MapTo_GPSSatellite(self, data): + prn = data['prn'] + doppler = dopplerOF.fromMapForm(data['doppler']) + amplitude = amplitudeOF.fromMapForm(data['amplitude']) + l1caEnabled = data['l1caEnabled'] + l2cEnabled = data['l2cEnabled'] + l1caMessage = messageOF.fromMapForm(data['l1caMessage']) + l2cMessage = messageOF.fromMapForm(data['l2cMessage']) + clCodeType = data['l2clCodeType'] + codeDopplerIgnored = data['codeDopplerIgnored'] + satellite = GPSSatellite(prn) + satellite.setAmplitude(amplitude) + satellite.setDoppler(doppler) + satellite.setL1CAEnabled(l1caEnabled) + satellite.setL2CEnabled(l2cEnabled) + satellite.setL1CAMessage(l1caMessage) + satellite.setL2CMessage(l2cMessage) + satellite.setL2CLCodeType(clCodeType) + satellite.setCodeDopplerIgnored(codeDopplerIgnored) + return satellite + +factoryObject = ObjectFactory() diff --git a/peregrine/iqgen/bits/satellite_gps.py b/peregrine/iqgen/bits/satellite_gps.py new file mode 100644 index 0000000..410804d --- /dev/null +++ b/peregrine/iqgen/bits/satellite_gps.py @@ -0,0 +1,227 @@ +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Valeri Atamaniouk +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + +""" +The :mod:`peregrine.iqgen.bits.satellite_gps` module contains classes and +functions related to GPS satellite configuration. + +""" +import peregrine.iqgen.bits.signals as signals +from peregrine.iqgen.bits.message_const import Message +from peregrine.iqgen.bits.prn_gps_l1ca import PrnCode as GPS_L1CA_Code +from peregrine.iqgen.bits.prn_gps_l2c import PrnCode as GPS_L2C_Code +from peregrine.iqgen.bits.satellite_base import Satellite + +import numpy + +DEFAULT_MESSAGE = Message(1) + + +class GPSSatellite(Satellite): + ''' + GPS satellite object. + ''' + + def __init__(self, prnNo): + ''' + Constructs satellite object + + Parameters + ---------- + prnNo : int + GPS satellite number for selecting PRN. + ''' + super(GPSSatellite, self).__init__("GPS{}".format(prnNo)) + self.prn = prnNo + self.l2clCodeType = '01' + self.l1caCode = GPS_L1CA_Code(prnNo) + self.l2cCode = GPS_L2C_Code(prnNo, self.l2clCodeType) + self.l1caEnabled = False + self.l2cEnabled = False + self.l1caMessage = DEFAULT_MESSAGE + self.l2cMessage = DEFAULT_MESSAGE + self.time0S = 0. + self.pr0M = 0. + self.phaseShift = 0. + + def setL1CAEnabled(self, enableFlag): + ''' + Enables or disable GPS L1 C/A sample generation + + Parameters + ---------- + enableFlag : boolean + Flag to enable (True) or disable (False) GPS L1 C/A samples + ''' + self.l1caEnabled = enableFlag + + def isL1CAEnabled(self): + ''' + Tests if L1 C/A signal generation is enabled + + Returns + ------- + bool + True, when L1 C/A signal generation is enabled, False otherwise + ''' + return self.l1caEnabled + + def setL2CEnabled(self, enableFlag): + ''' + Enables or disable GPS L2 C sample generation + + Parameters + ---------- + enableFlag : boolean + Flag to enable (True) or disable (False) GPS L2 C samples + ''' + self.l2cEnabled = enableFlag + + def isL2CEnabled(self): + ''' + Tests if L2 C signal generation is enabled + + Returns + ------- + bool + True, when L2 C signal generation is enabled, False otherwise + ''' + return self.l2cEnabled + + def setL2CLCodeType(self, clCodeType): + if self.l2clCodeType != clCodeType: + self.l2cCode = GPS_L2C_Code(self.prn, clCodeType) + self.l2clCodeType = clCodeType + + def getL2CLCodeType(self): + return self.l2clCodeType + + def setL1CAMessage(self, message): + ''' + Configures data source for L1 C/A signal. + + Parameters + ---------- + message : object + Message object that will provide symbols for L1 C/A signal. + ''' + self.l1caMessage = message + + def setL2CMessage(self, message): + ''' + Configures data source for L2 C signal. + + Parameters + ---------- + message : object + Message object that will provide symbols for L2 C signal. + ''' + self.l2cMessage = message + + def getL1CAMessage(self): + ''' + Returns configured message object for GPS L1 C/A band + + Returns + ------- + object + L1 C/A message object + ''' + return self.l1caMessage + + def getL2CMessage(self): + ''' + Returns configured message object for GPS L2 C band + + Returns + ------- + object + L2 C message object + ''' + return self.l2cMessage + + def getBatchSignals(self, userTimeAll_s, samples, outputConfig, debug): + ''' + Generates signal samples. + + Parameters + ---------- + userTimeAll_s : numpy.ndarray(n_samples, dtype=numpy.float64) + Vector of observer's timestamps in seconds for the interval start. + samples : numpy.ndarray((4, n_samples)) + Array to which samples are added. + outputConfig : object + Output configuration object. + debug : bool + Debug flag + + Returns + ------- + list + Debug information + ''' + result = [] + if (self.l1caEnabled): + intermediateFrequency_hz = outputConfig.GPS.L1.INTERMEDIATE_FREQUENCY_HZ + frequencyIndex = outputConfig.GPS.L1.INDEX + values = self.doppler.computeBatch(userTimeAll_s, + self.amplitude, + signals.GPS.L1CA, + intermediateFrequency_hz, + self.l1caMessage, + self.l1caCode, + outputConfig, + debug) + numpy.add(samples[frequencyIndex], + values[0], + out=samples[frequencyIndex]) + debugData = {'type': "GPSL1", 'doppler': values[1]} + result.append(debugData) + if (self.l2cEnabled): + intermediateFrequency_hz = outputConfig.GPS.L2.INTERMEDIATE_FREQUENCY_HZ + frequencyIndex = outputConfig.GPS.L2.INDEX + values = self.doppler.computeBatch(userTimeAll_s, + self.amplitude, + signals.GPS.L2C, + intermediateFrequency_hz, + self.l2cMessage, + self.l2cCode, + outputConfig, + debug) + numpy.add(samples[frequencyIndex], + values[0], + out=samples[frequencyIndex]) + debugData = {'type': "GPSL2", 'doppler': values[1]} + result.append(debugData) + return result + + def isBandEnabled(self, bandIndex, outputConfig): + ''' + Checks if particular band is supported and enabled. + + Parameters + ---------- + bandIndex : int + Signal band index + outputConfig : object + Output configuration + + Returns: + bool + True, if the band is supported and enabled; False otherwise. + ''' + result = None + if bandIndex == outputConfig.GPS.L1.INDEX: + result = self.isL1CAEnabled() + elif bandIndex == outputConfig.GPS.L2.INDEX: + result = self.isL2CEnabled() + else: + result = False + return result diff --git a/peregrine/iqgen/bits/signals.py b/peregrine/iqgen/bits/signals.py new file mode 100644 index 0000000..246fd76 --- /dev/null +++ b/peregrine/iqgen/bits/signals.py @@ -0,0 +1,171 @@ +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Valeri Atamaniouk +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + +""" +The :mod:`peregrine.iqgen.bits.signals` module contains classes and functions +related to common satellite signal definitions and operations. + +""" + +import scipy.constants + + +class GPS: + ''' + GPS signal parameters and utilities. + ''' + + @staticmethod + def _calcDopplerShiftHz(frequency_hz, distance_m, velocity_mps): + ''' + Utility to compute doppler shift from ditance and velocity for a band + frequency. + + Parameters + ---------- + frequency_hz : float + Band frequency in hertz + distance_m : float + Distance to satellite in meters + velocity_m : float + Satellite velocity in meters per second. + + Return + ------ + float + Doppler shift in hertz + ''' + doppler_hz = -velocity_mps * frequency_hz / scipy.constants.c + return doppler_hz + + class L1CA: + ''' + GPS L1 C/A parameters and utilities. + ''' + SYMBOL_RATE_HZ = 50 + CENTER_FREQUENCY_HZ = 1575.42e6 + CODE_CHIP_RATE_HZ = 1023000 + CHIP_TO_SYMBOL_DIVIDER = 20460 + + @staticmethod + def calcDopplerShiftHz(distance_m, velocity_mps): + ''' + Converts relative speed into doppler value for GPS L1 C/A band. + + Parameters + ---------- + distance_m : float + Distance in meters + velocity_mps : float + Relative speed in meters per second. + + Returns + ------- + float + Doppler shift in Hz. + ''' + return GPS._calcDopplerShiftHz(GPS.L1CA.CENTER_FREQUENCY_HZ, distance_m, velocity_mps) + + @staticmethod + def getSymbolIndex(svTime_s): + ''' + Computes symbol index. + + Parameters + ---------- + svTime_s : float + SV time in seconds + + Returns + ------- + long + Symbol index + ''' + return long(svTime_s * GPS.L1CA.SYMBOL_RATE_HZ) + + @staticmethod + def getCodeChipIndex(svTime_s): + ''' + Computes code chip index. + + Parameters + ---------- + svTime_s : float + SV time in seconds + + Returns + ------- + long + Code chip index + ''' + return long(svTime_s * GPS.L1CA.CODE_CHIP_RATE_HZ) + + class L2C: + ''' + GPS L2 C parameters and utilities. + ''' + + SYMBOL_RATE_HZ = 50 + CENTER_FREQUENCY_HZ = 1227.60e6 + CODE_CHIP_RATE_HZ = 1023000 + CHIP_TO_SYMBOL_DIVIDER = 20460 + + @staticmethod + def calcDopplerShiftHz(distance_m, velocity_mps): + ''' + Converts relative speed into doppler value for GPS L2 C band. + + Parameters + ---------- + distance_m : float + Distance in meters + velocity_mps : float + Relative speed in meters per second. + + Returns + ------- + float + Doppler shift in Hz. + ''' + return GPS._calcDopplerShiftHz(GPS.L2C.CENTER_FREQUENCY_HZ, distance_m, velocity_mps) + + @staticmethod + def getSymbolIndex(svTime_s): + ''' + Computes symbol index. + + Parameters + ---------- + svTime_s : float + SV time in seconds + + Returns + ------- + long + Symbol index + ''' + return long(svTime_s * GPS.L2C.SYMBOL_RATE_HZ) + + @staticmethod + def getCodeChipIndex(svTime_s): + ''' + Computes code chip index. + + Parameters + ---------- + svTime_s : float + SV time in seconds + + Returns + ------- + long + Code chip index + ''' + return long(svTime_s * GPS.L2C.CODE_CHIP_RATE_HZ) diff --git a/peregrine/iqgen/bits/tcxo_base.py b/peregrine/iqgen/bits/tcxo_base.py new file mode 100644 index 0000000..3bd0585 --- /dev/null +++ b/peregrine/iqgen/bits/tcxo_base.py @@ -0,0 +1,46 @@ +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Valeri Atamaniouk +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + +""" +The :mod:`peregrine.iqgen.bits.tcxo_base` module contains base class definitions +for TCXO control. + +""" + + +class TCXOBase(object): + ''' + Base class for TCXO control. The class computes time shifts of TCXO depending + on external conditions like temperature, vibration, etc. + ''' + + def __init__(self): + super(TCXOBase, self).__init__() + + def computeTcxoTime(self, fromSample, toSample, outputConfig): + ''' + Method generates time vector for the given sample index range depending on + TCXO behaviour. + + Parameters + ---------- + fromSample : int + Index of the first sample. + toSample: int + Index of the last sample plus 1. + outputConfig : object + Output configuration + + Returns + ------- + numpy.ndarray(shape=(toSample - fromSample), dtype=numpy.float) + Vector of the shifted time stamps for the given TCXO controller. + ''' + raise NotImplementedError() diff --git a/peregrine/iqgen/bits/tcxo_factory.py b/peregrine/iqgen/bits/tcxo_factory.py new file mode 100644 index 0000000..6501654 --- /dev/null +++ b/peregrine/iqgen/bits/tcxo_factory.py @@ -0,0 +1,68 @@ +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Valeri Atamaniouk +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + +""" +The :mod:`peregrine.iqgen.bits.amplitude_factory` module contains classes and +functions related to object factory for amplitude objects. + +""" + +from peregrine.iqgen.bits.tcxo_poly import TCXOPoly as PolyTcxo +from peregrine.iqgen.bits.tcxo_sine import TCXOSine as SineTcxo + + +class ObjectFactory(object): + ''' + Object factory for amplitude objects. + ''' + + def __init__(self): + super(ObjectFactory, self).__init__() + + def toMapForm(self, obj): + t = type(obj) + if t is PolyTcxo: + return self.__PolyTcxo_ToMap(obj) + elif t is SineTcxo: + return self.__SineTcxo_ToMap(obj) + else: + raise ValueError("Invalid object type") + + def fromMapForm(self, data): + t = data['type'] + if t == 'PolyTcxo': + return self.__MapTo_PolyTcxo(data) + elif t == 'SineTcxo': + return self.__MapTo_SineTcxo(data) + else: + raise ValueError("Invalid object type") + + def __PolyTcxo_ToMap(self, obj): + data = {'type': 'PolyTcxo', 'coeffs': obj.coeffs} + return data + + def __SineTcxo_ToMap(self, obj): + data = {'type': 'SineTcxo', + 'initial_ppm': obj.initial_ppm, + 'amplitude_ppm': obj.amplitude_ppm, + 'period_s': obj.period_s} + return data + + def __MapTo_PolyTcxo(self, data): + coeffs = data['coeffs'] + return PolyTcxo(coeffs) + + def __MapTo_SineTcxo(self, data): + initial_ppm = data['initial_ppm'] + amplitude_ppm = data['amplitude_ppm'] + period_s = data['period_s'] + return SineTcxo(initial_ppm, amplitude_ppm, period_s) + +factoryObject = ObjectFactory() diff --git a/peregrine/iqgen/bits/tcxo_poly.py b/peregrine/iqgen/bits/tcxo_poly.py new file mode 100644 index 0000000..f7fcbe5 --- /dev/null +++ b/peregrine/iqgen/bits/tcxo_poly.py @@ -0,0 +1,95 @@ +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Valeri Atamaniouk +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + +""" +The :mod:`peregrine.iqgen.bits.tcxo_poly` module contains class definitions +for TCXO control that can describe TCXO frequency drift as a polynome. + +""" + +from peregrine.iqgen.bits.tcxo_base import TCXOBase +import numpy + + +class TCXOPoly(TCXOBase): + ''' + Polynomial TCXO control class. + ''' + + def __init__(self, coeffs): + ''' + Constructs TCXO control object. + + Parameters + ---------- + coeffs : array-like + Coefficients for TCXO polynome. These coeffificens define a TCXO drift + over time in ppm. + ''' + super(TCXOPoly, self).__init__() + self.coeffs = coeffs[:] + if coeffs: + # Recompute drift coefficients from speed of drift into distance of drift + new_coeffs = [] + power_c = len(coeffs) + for idx, val in enumerate(coeffs): + power = power_c - idx + new_coeffs.append(val * 1e-6 / power) + new_coeffs.append(0) + self.poly = numpy.poly1d(new_coeffs) + else: + self.poly = None + + def __str__(self, *args, **kwargs): + ''' + Provides string representation of the object + ''' + return "TCXOPoly: coeffs=%s" % str(self.coeffs) + + def __repr__(self): + ''' + Provides string representation of the object + ''' + return "TCXOPoly(%s)" % repr(self.coeffs) + + def computeTcxoTime(self, fromSample, toSample, outputConfig): + ''' + Method generates time vector for the given sample index range depending on + TCXO behaviour. + + Parameters + ---------- + fromSample : int + Index of the first sample. + toSample: int + Index of the last sample plus 1. + outputConfig : object + Output configuration + + Returns + ------- + numpy.ndarray(shape=(toSample - fromSample), dtype=numpy.float) + Vector of the shifted time stamps for the given TCXO controller. + ''' + poly = self.poly + + if poly: + time0_s = fromSample / outputConfig.SAMPLE_RATE_HZ + timeX_s = toSample / outputConfig.SAMPLE_RATE_HZ + timeAll_s = numpy.linspace(time0_s, + timeX_s, + toSample - fromSample, + endpoint=False, + dtype=numpy.float) + result = poly(timeAll_s) + else: + result = None + + return result diff --git a/peregrine/iqgen/bits/tcxo_sine.py b/peregrine/iqgen/bits/tcxo_sine.py new file mode 100644 index 0000000..cc0b3b3 --- /dev/null +++ b/peregrine/iqgen/bits/tcxo_sine.py @@ -0,0 +1,101 @@ +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Valeri Atamaniouk +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + +""" +The :mod:`peregrine.iqgen.bits.tcxo_sine` module contains class definitions +for TCXO control that can describe TCXO frequency drift as a periodic (sine) +function. + +""" + +from peregrine.iqgen.bits.tcxo_base import TCXOBase +import numpy +import scipy.constants + + +class TCXOSine(TCXOBase): + ''' + Sine TCXO control class. + ''' + + def __init__(self, initial_ppm, amplitude_ppm, period_s): + ''' + Constructs TCXO control object. + + Parameters + ---------- + initial_ppm : float + Initial drift in ppm + amplitude_ppm : float + Drift amplitude in ppm + period_s : float + Drift period in seconds + ''' + super(TCXOSine, self).__init__() + + self.initial_ppm = initial_ppm + self.amplitude_ppm = amplitude_ppm + self.period_s = period_s + self.c0 = -2. * scipy.constants.pi * amplitude_ppm * 1e-6 + self.c1 = 2. * scipy.constants.pi / period_s + self.c2 = initial_ppm * 1e-6 + + def __str__(self, *args, **kwargs): + ''' + Provides string representation of the object + ''' + return "TCXOSine: initial_ppm=%f amplitude_ppm=%f period_s=%f" % \ + (self.initial_ppm, self.amplitude_ppm, self.period_s) + + def __repr__(self): + ''' + Provides string representation of the object + ''' + return "TCXOSine(%f, %f, %f)" % \ + (self.initial_ppm, self.amplitude_ppm, self.period_s) + + def computeTcxoTime(self, fromSample, toSample, outputConfig): + ''' + Method generates time vector for the given sample index range depending on + TCXO behaviour. + + Parameters + ---------- + fromSample : int + Index of the first sample. + toSample: int + Index of the last sample plus 1. + outputConfig : object + Output configuration + + Returns + ------- + numpy.ndarray(shape=(toSample - fromSample), dtype=numpy.float) + Vector of the shifted time stamps for the given TCXO controller. + ''' + c0 = self.c0 + c1 = self.c1 + c2 = self.c2 + time0_s = fromSample / outputConfig.SAMPLE_RATE_HZ + timeX_s = toSample / outputConfig.SAMPLE_RATE_HZ + + timeAll_s = numpy.linspace(time0_s * c1, + timeX_s * c1, + toSample - fromSample, + endpoint=False, + dtype=numpy.float) + + result = numpy.cos(timeAll_s) + result += -1. + result *= c0 + if c2: + result += timeAll_s * c2 + + return result diff --git a/peregrine/iqgen/generate.py b/peregrine/iqgen/generate.py new file mode 100644 index 0000000..4d9a826 --- /dev/null +++ b/peregrine/iqgen/generate.py @@ -0,0 +1,598 @@ +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Valeri Atamaniouk +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + +import sys +import traceback + +""" +The :mod:`peregrine.iqgen.generate` module contains classes and functions +related to main loop of samples generation. + +""" +from peregrine.iqgen.bits.filter_lowpass import LowPassFilter +from peregrine.iqgen.bits.filter_bandpass import BandPassFilter + +from peregrine.iqgen.bits import signals + +import logging +import scipy +import numpy +import time + +import multiprocessing + +logger = logging.getLogger(__name__) + + +class Task(object): + ''' + Period computation task. This object performs a batch computation of signal + in the specified range. + ''' + + def __init__(self, + outputConfig, + signalSources, + noiseSigma, + tcxo, + signalFilters, + generateDebug): + ''' + Parameters + ---------- + outputConfig : object + Output profile + signalSources : array-like + List of satellites + noiseSigma : float + Noise sigma value + tcxo : object + TCXO control object + signalFilters : array-like + Output signal filter objects + generateDebug : bool + Flag if additional debug output is required + ''' + + self.outputConfig = outputConfig + self.signalSources = signalSources + self.signalFilters = signalFilters + self.generateDebug = generateDebug + self.noiseSigma = noiseSigma + self.tcxo = tcxo + self.signals = scipy.ndarray(shape=(4, outputConfig.SAMPLE_BATCH_SIZE), + dtype=numpy.float) + self.noise = self.createNoise(outputConfig.SAMPLE_BATCH_SIZE) + self.nSamples = outputConfig.SAMPLE_BATCH_SIZE + + def update(self, userTime0_s, nSamples, firstSampleIndex): + ''' + Configure object for the next batch generation + + The method stores parameters for the generation step and updates internal + arrays to match output shape. + + Parameters + ---------- + userTime0_s : float + Time of the interval start in seconds + nSamples : long + Number of samples in the interval + firstSampleIndex : long + Index of the first sample + ''' + self.userTime0_s = userTime0_s + self.firstSampleIndex = firstSampleIndex + + if (self.nSamples != nSamples): + newSignals = numpy.ndarray((4, nSamples), dtype=float) + newNoise = self.createNoise(nSamples) + self.nSamples = nSamples + self.signals = newSignals + self.noise = newNoise + + def createNoise(self, nSamples): + ''' + Generate noise array for a given noise sigma. + + Parameters + ---------- + nSamples : int + Length of the samples vectors + + Returns + ------- + numpy.ndarray(shape=(4, nSamples), dtype=numpy.float) + Noise values + ''' + noiseSigma = self.noiseSigma + if noiseSigma is not None: + # Initialize signal array with noise + noiseType = 1 + if noiseType == 1: + noise = noiseSigma * scipy.randn(4, nSamples) + else: + noise = numpy.random.normal(loc=0., + scale=noiseSigma, + size=(4, nSamples)) + # print self.noise + else: + noise = None + return noise + + def perform(self): + outputConfig = self.outputConfig + signalSources = self.signalSources + signalFilters = self.signalFilters + tcxo = self.tcxo + firstSampleIndex = self.firstSampleIndex + finalSampleIndex = firstSampleIndex + self.nSamples + + generateDebug = self.generateDebug + + userTime0_s = self.userTime0_s + userTimeX_s = userTime0_s + float(self.nSamples) / \ + float(outputConfig.SAMPLE_RATE_HZ) + userTimeAll_s = scipy.linspace(userTime0_s, + userTimeX_s, + self.nSamples, + endpoint=False) + + if tcxo: + tcxoTimeDrift_s = tcxo.computeTcxoTime(firstSampleIndex, + finalSampleIndex, + outputConfig) + if tcxoTimeDrift_s: + userTimeAll_s += tcxoTimeDrift_s + + noise = self.noise + sigs = self.signals + sigs.fill(0.) + if noise is not None: + # Initialize signal array with noise + sigs += noise + + if generateDebug: + signalData = [] + debugData = {'time': userTimeAll_s, 'signalData': signalData} + else: + debugData = None + + # Sum up signals for all SVs + for signalSource in signalSources: + # Add signal from source (satellite) to signal accumulator + t = signalSource.getBatchSignals(userTimeAll_s, + sigs, + outputConfig, + generateDebug) + # Debugging output + if generateDebug: + svDebug = {'name': signalSource.getSvName(), 'data': t} + signalData.append(svDebug) + t = None + + if signalFilters is list: + # Filter signal values through LPF, BPF or another + for i in range(len(self.filters)): + filterObject = signalFilters[i] + if filterObject is not None: + sigs[i][:] = filterObject.filter(sigs[i]) + + inputParams = (self.userTime0_s, self.nSamples, self.firstSampleIndex) + return (inputParams, sigs, debugData) + + +class Worker(multiprocessing.Process): + + def __init__(self, + outputConfig, + signalSources, + noiseSigma, + tcxo, + signalFilters, + generateDebug): + super(Worker, self).__init__() + self.queueIn = multiprocessing.Queue() + self.queueOut = multiprocessing.Queue() + self.totalWaitTime_s = 0. + self.totalExecTime_s = 0. + self.outputConfig = outputConfig + self.signalSources = signalSources + self.noiseSigma = noiseSigma + self.tcxo = tcxo + self.signalFilters = signalFilters + self.generateDebug = generateDebug + + def run(self): + task = Task(self.outputConfig, + self.signalSources, + noiseSigma=self.noiseSigma, + tcxo=self.tcxo, + signalFilters=self.signalFilters, + generateDebug=self.generateDebug) + + while True: + opStartTime_s = time.clock() + inputRequest = self.queueIn.get() + if inputRequest is None: + # EOF reached + break + (userTime0_s, nSamples, firstSampleIndex) = inputRequest + # print "Received params", userTime0_s, nSamples, firstSampleIndex + opDuration_s = time.clock() - opStartTime_s + self.totalWaitTime_s += opDuration_s + startTime_s = time.clock() + try: + task.update(userTime0_s, nSamples, firstSampleIndex) + result = task.perform() + import copy + result = copy.deepcopy(result) + self.queueOut.put(result) + except: + exType, exValue, exTraceback = sys.exc_info() + traceback.print_exception( + exType, exValue, exTraceback, file=sys.stderr) + self.queueOut.put(None) + self.queueIn.close() + self.queueOut.close() + sys.exit(1) + duration_s = time.clock() - startTime_s + self.totalExecTime_s += duration_s + statistics = (self.totalWaitTime_s, self.totalExecTime_s) + self.queueOut.put(statistics) + self.queueIn.close() + self.queueOut.close() + sys.exit(0) + + +def computeTimeIntervalS(outputConfig): + ''' + Helper for computing generation interval duration in seconds. + + Parameters + ---------- + outputConfig : object + Output configuration. + + Returns + ------- + float + Generation interval duration in seconds + ''' + deltaTime_s = float(outputConfig.SAMPLE_BATCH_SIZE) / \ + outputConfig.SAMPLE_RATE_HZ + return deltaTime_s + + +def generateSamples(outputFile, + sv_list, + encoder, + time0S, + nSamples, + outputConfig, + SNR=None, + tcxo=None, + filterType="none", + logFile=None, + threadCount=0, + pbar=None): + ''' + Generates samples. + + Parameters + ---------- + fileName : string + Output file name. + sv_list : list + List of configured satellite objects. + encoder : Encoder + Output encoder object. + time0S : float + Time epoch for the first sample. + nSamples : long + Total number of samples to generate. + outputConfig : object + Output parameters + SNR : float, optional + When specified, adds random noise to the output. + tcxo : object, optional + When specified, controls TCXO drift + filterType : string, optional + Controls IIR/FIR signal post-processing. Disabled by default. + debugLog : bool, optional + Control generation of additional debug output. Disabled by default. + ''' + + # + # Print out parameters + # + print "Generating samples, sample rate={} Hz, interval={} seconds, SNR={}".format( + outputConfig.SAMPLE_RATE_HZ, nSamples / outputConfig.SAMPLE_RATE_HZ, SNR) + print "Jobs: ", threadCount + + _t0 = time.clock() + _count = 0l + + # Check which bands are enabled, configure band-specific parameters + bands = [outputConfig.GPS.L1, outputConfig.GPS.L2] # Supported bands + lpf = [None] * len(bands) + bandsEnabled = [False] * len(bands) + + bandPass = False + lowPass = False + if filterType == 'lowpass': + lowPass = True + elif filterType == 'bandpass': + bandPass = True + elif filterType == 'none': + pass + else: + raise ValueError("Invalid filter type %s" % repr(filter)) + + for band in bands: + for sv in sv_list: + bandsEnabled[band.INDEX] |= sv.isBandEnabled(band.INDEX, outputConfig) + + filterObject = None + if lowPass: + filterObject = LowPassFilter(outputConfig, + band.INTERMEDIATE_FREQUENCY_HZ) + elif bandPass: + filterObject = BandPassFilter(outputConfig, + band.INTERMEDIATE_FREQUENCY_HZ) + if filterObject: + lpf[band.INDEX] = filterObject + logger.debug("Band %d filter NBW is %s" % + (band.INDEX, str(filterObject))) + + if SNR is not None: + sourcePower = 0. + for sv in sv_list: + svMeanPower = sv.getAmplitude().computeMeanPower() + sourcePower += svMeanPower + logger.debug("[%s] Estimated mean power is %f" % + (sv.getSvName(), svMeanPower)) + meanPower = sourcePower / len(sv_list) + meanAmplitude = scipy.sqrt(meanPower) + logger.debug("Estimated total signal power is %f, mean %f, mean amplitude %f" % + (sourcePower, meanPower, meanAmplitude)) + + # Nsigma and while noise amplitude computation: check if the Nsigma is + # actually a correct value for white noise with normal distribution. + + # Number of samples for 1023 MHz + freqTimesTau = outputConfig.SAMPLE_RATE_HZ / 1.023e6 + noiseVariance = freqTimesTau * meanPower / (4. * 10. ** (float(SNR) / 10.)) + noiseSigma = numpy.sqrt(noiseVariance) + logger.info("Selected noise sigma %f (variance %f) for SNR %f" % + (noiseSigma, noiseVariance, float(SNR))) + + else: + noiseVariance = None + noiseSigma = None + logger.info("SNR is not provided, noise is not generated.") + + # + # Print out SV parameters + # + for _sv in sv_list: + _svNo = _sv.getSvName() + _amp = _sv.amplitude + _svTime0_s = 0 + _dist0_m = _sv.doppler.computeDistanceM(_svTime0_s) + _speed_mps = _sv.doppler.computeSpeedMps(_svTime0_s) + _bit = signals.GPS.L1CA.getSymbolIndex(_svTime0_s) + _c1 = signals.GPS.L1CA.getCodeChipIndex(_svTime0_s) + _c2 = signals.GPS.L2C.getCodeChipIndex(_svTime0_s) + _d1 = signals.GPS.L1CA.calcDopplerShiftHz(_dist0_m, _speed_mps) + _d2 = signals.GPS.L2C.calcDopplerShiftHz(_dist0_m, _speed_mps) + svMeanPower = _sv.getAmplitude().computeMeanPower() + # SNR for a satellite. Depends on sampling rate. + if noiseVariance: + svSNR = svMeanPower / (4. * noiseVariance) * freqTimesTau + else: + svSNR = 1e6 + svSNR_db = 10. * numpy.log10(svSNR) + # Filters lower the power according to their attenuation levels + l1FA_db = lpf[0].getPassBandAtt() if lpf[0] else 0. + l2FA_db = lpf[1].getPassBandAtt() if lpf[1] else 0. + # CNo for L1 + svCNoL1 = svSNR_db + 10. * numpy.log10(1.023e6) - l1FA_db + # CNo for L2, half power used (-3dB) + svCNoL2 = svSNR_db + 10. * numpy.log10(1.023e6) - 3. - l2FA_db + + print "{} = {{".format(_svNo) + print " .amplitude: {}".format(_amp) + print " .doppler: {}".format(_sv.doppler) + print " .l1_message: {}".format(_sv.getL1CAMessage()) + print " .l2_message: {}".format(_sv.getL2CMessage()) + print " .l2_cl_type: {}".format(_sv.getL2CLCodeType()) + print " .SNR (dBHz): {}".format(svSNR_db) + print " .L1 CNo: {}".format(svCNoL1) + print " .L2 CNo: {}".format(svCNoL2) + print " .epoc:" + print " .distance: {} m".format(_dist0_m) + print " .speed: {} m/s".format(_speed_mps) + print " .symbol: {}".format(_bit) + print " .l1_doppler: {} hz".format(_d1) + print " .l2_doppler: {} hz".format(_d2) + print " .l1_chip: {}".format(_c1) + print " .l2_chip: {}".format(_c2) + print "}" + + userTime_s = float(time0S) + + deltaUserTime_s = computeTimeIntervalS(outputConfig) + debugFlag = logFile is not None + + if debugFlag: + logFile.write("Index,Time") + for sv in sv_list: + svName = sv.getSvName() + if sv.isL1CAEnabled(): + logFile.write(",%s/L1/doppler" % svName) + if sv.isL2CEnabled(): + logFile.write(",%s/L2/doppler" % svName) + # End of line + logFile.write("\n") + + if threadCount > 0: + workerPool = [Worker(outputConfig, + sv_list, + noiseSigma, + tcxo, + lpf, + debugFlag) for _ in range(threadCount)] + + for worker in workerPool: + worker.start() + maxTaskListSize = threadCount * 2 + else: + workerPool = None + task = Task(outputConfig, + sv_list, + noiseSigma=noiseSigma, + tcxo=tcxo, + signalFilters=lpf, + generateDebug=debugFlag) + maxTaskListSize = 1 + + workerPutIndex = 0 + workerGetIndex = 0 + activeTasks = 0 + + totalSampleCounter = 0 + taskQueuedCounter = 0 + taskReceivedCounter = 0 + + totalEncodeTime_s = 0. + totalWaitTime_s = 0. + + while True: + while activeTasks < maxTaskListSize and totalSampleCounter < nSamples: + # We have space in the task backlog and not all batchIntervals are issued + userTime0_s = userTime_s + userTimeX_s = userTime_s + deltaUserTime_s + sampleCount = outputConfig.SAMPLE_BATCH_SIZE + + if totalSampleCounter + sampleCount > nSamples: + # Last interval may contain less than full batch size of samples + sampleCount = nSamples - totalSampleCounter + userTimeX_s = userTime0_s + float(sampleCount) / \ + outputConfig.SAMPLE_RATE_HZ + + params = (userTime0_s, sampleCount, totalSampleCounter) + # print ">>> ", userTime0_s, sampleCount, totalSampleCounter, + # workerPutIndex + if workerPool is not None: + workerPool[workerPutIndex].queueIn.put(params) + workerPutIndex = (workerPutIndex + 1) % threadCount + else: + task.update(userTime0_s, sampleCount, totalSampleCounter) + activeTasks += 1 + + # Update parameters for the next batch interval + userTime_s = userTimeX_s + totalSampleCounter += sampleCount + taskQueuedCounter += 1 + + # What for the data only if we have something to wait + if taskReceivedCounter == taskQueuedCounter and \ + totalSampleCounter == nSamples: + # No more tasks to issue to generator + # No more tasks to wait + break + + try: + if workerPool is not None: + # Wait for the first task + worker = workerPool[workerGetIndex] + waitStartTime_s = time.time() + # print "waiting data from worker", workerGetIndex + result = worker.queueOut.get() + # print "Data received from worker", workerGetIndex + workerGetIndex = (workerGetIndex + 1) % threadCount + waitDuration_s = time.time() - waitStartTime_s + totalWaitTime_s += waitDuration_s + else: + result = task.perform() + except: + exType, exValue, exTraceback = sys.exc_info() + traceback.print_exception(exType, exValue, exTraceback, file=sys.stderr) + result = None + taskReceivedCounter += 1 + activeTasks -= 1 + + if result is None: + print "Error in processor; aborting." + break + + (inputParams, signalSamples, debugData) = result + (_userTime0_s, _sampleCount, _firstSampleIndex) = inputParams + # print "<<< ", _userTime0_s, _sampleCount, _firstSampleIndex + + if logFile is not None: + # Data from all satellites is collected. Now we can dump the debug matrix + + userTimeAll_s = debugData['time'] + signalData = debugData['signalData'] + for smpl_idx in range(_sampleCount): + logFile.write("{},{}".format(_firstSampleIndex + smpl_idx, + userTimeAll_s[smpl_idx])) + for svIdx in range(len(signalData)): + # signalSourceName = signalData[svIdx]['name'] + signalSourceData = signalData[svIdx]['data'] + for band in signalSourceData: + # bandType = band['type'] + doppler = band['doppler'] + logFile.write(",{}".format(doppler[smpl_idx])) + # End of line + logFile.write("\n") + + encodeStartTime_s = time.time() + # Feed data into encoder + encodedSamples = encoder.addSamples(signalSamples) + signalSamples = None + + if len(encodedSamples) > 0: + _count += len(encodedSamples) + encodedSamples.tofile(outputFile) + encodedSamples = None + encodeDuration_s = time.time() - encodeStartTime_s + totalEncodeTime_s += encodeDuration_s + + if pbar: + pbar.update(_firstSampleIndex + _sampleCount) + + logger.debug("MAIN: Encode duration: %f" % totalEncodeTime_s) + logger.debug("MAIN: wait duration: %f" % totalWaitTime_s) + + encodedSamples = encoder.flush() + if len(encodedSamples) > 0: + encodedSamples.tofile(outputFile) + + if debugFlag: + logFile.close() + + if workerPool is not None: + for worker in workerPool: + worker.queueIn.put(None) + for worker in workerPool: + try: + statistics = worker.queueOut.get(timeout=2) + print "Statistics:", statistics + except: + exType, exValue, exTraceback = sys.exc_info() + traceback.print_exception( + exType, exValue, exTraceback, file=sys.stderr) + worker.queueIn.close() + worker.queueOut.close() + worker.terminate() + worker.join() diff --git a/peregrine/iqgen/if_iface.py b/peregrine/iqgen/if_iface.py new file mode 100644 index 0000000..b432b35 --- /dev/null +++ b/peregrine/iqgen/if_iface.py @@ -0,0 +1,222 @@ +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Valeri Atamaniouk +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + +""" +The :mod:`peregrine.iqgen.if_iface` module contains classes and functions +related to radio interface parameters + +""" + +from peregrine.defaults import freq_profile_peregrine + + +class LowRateConfig(object): + ''' + Output control configuration for quick tests. + + Attributes + ---------- + NAME : string + Configuration name + SAMPLE_RATE_HZ : float + Sample rate in hertz for data generation. + SAMPLE_BATCH_SIZE : int + Size of the sample batch in samples. + GPS : object + GPS band information + Galileo : object + Galileo band information + Beidou : object + Beidou band information + Glonass : object + Glonass band information + ''' + NAME = "Low rate configuration for fast tests" + SAMPLE_RATE_HZ = 24.84375e5 + SAMPLE_BATCH_SIZE = 100000 + + class GPS(object): + + class L1(object): + INTERMEDIATE_FREQUENCY_HZ = 14.58e5 + INDEX = 0 + + class L2(object): + INTERMEDIATE_FREQUENCY_HZ = 7.4e+5 + INDEX = 1 + + class Glonass(object): + + class L1(object): + INTERMEDIATE_FREQUENCY_HZ = 12e5 + INDEX = 1 + + class L2(object): + INTERMEDIATE_FREQUENCY_HZ = 11e5 + INDEX = 2 + + class Galileo(object): + + class E1(object): + INTERMEDIATE_FREQUENCY_HZ = 14.58e5 + INDEX = 0 + + class E6(object): + INTERMEDIATE_FREQUENCY_HZ = 43.75e5 + INDEX = 2 + + class E5b(object): + INTERMEDIATE_FREQUENCY_HZ = 27.86e5 + INDEX = 3 + + class Beidou(object): + + class B1(object): + INTERMEDIATE_FREQUENCY_HZ = 28.902e5 + INDEX = 0 + + class B2: + INTERMEDIATE_FREQUENCY_HZ = 27.86e5 + INDEX = 3 + + class B3(object): + INTERMEDIATE_FREQUENCY_HZ = 33.52e5 + INDEX = 2 + + +class NormalRateConfig(object): + ''' + Output control configuration for normal tests. + + Attributes + ---------- + NAME : string + Configuration name + SAMPLE_RATE_HZ : float + Sample rate in hertz for data generation. + SAMPLE_BATCH_SIZE : int + Size of the sample batch in samples. + GPS : object + GPS band information + Galileo : object + Galileo band information + Beidou : object + Beidou band information + Glonass : object + Glonass band information + ''' + NAME = "Normal rate configuration equivalent to decimated data output" + SAMPLE_RATE_HZ = 24.84375e6 + SAMPLE_BATCH_SIZE = 100000 + + class GPS(object): + ''' + Parameters for GPS bands data generation. + ''' + class L1(object): + INTERMEDIATE_FREQUENCY_HZ = 14.58e+6 + INDEX = 0 + + class L2(object): + INTERMEDIATE_FREQUENCY_HZ = 7.4e+6 + INDEX = 1 + + class Glonass(object): + + class L1(object): + INTERMEDIATE_FREQUENCY_HZ = 12e6 + INDEX = 1 + + class L2(object): + INTERMEDIATE_FREQUENCY_HZ = 11e6 + INDEX = 2 + + class Galileo(object): + + class E1(object): + INTERMEDIATE_FREQUENCY_HZ = 14.58e6 + INDEX = 0 + + class E6(object): + INTERMEDIATE_FREQUENCY_HZ = 43.75e6 + INDEX = 2 + + class E5b(object): + INTERMEDIATE_FREQUENCY_HZ = 27.86e6 + INDEX = 3 + + class Beidou(object): + + class B1(object): + INTERMEDIATE_FREQUENCY_HZ = 28.902e6 + INDEX = 0 + + class B2: + INTERMEDIATE_FREQUENCY_HZ = 27.86e6 + INDEX = 3 + + class B3(object): + INTERMEDIATE_FREQUENCY_HZ = 33.52e6 + INDEX = 2 + + +class HighRateConfig(object): + ''' + Output control configuration for high data rate tests. + + Attributes + ---------- + NAME : string + Configuration name + SAMPLE_RATE_HZ : float + Sample rate in hertz for data generation. + SAMPLE_BATCH_SIZE : int + Size of the sample batch in samples. + GPS : object + GPS band information + ''' + NAME = "High rate configuration equivalent to full rate data output" + SAMPLE_RATE_HZ = 99.375e6 + SAMPLE_BATCH_SIZE = 100000 + + GPS = NormalRateConfig.GPS + Glonass = NormalRateConfig.Glonass + Galileo = NormalRateConfig.Galileo + Beidou = NormalRateConfig.Beidou + + +class CustomRateConfig(object): + ''' + Output control configuration for comparison tests. + + Attributes + ---------- + NAME : string + Configuration name + SAMPLE_RATE_HZ : float + Sample rate in hertz for data generation. + SAMPLE_BATCH_SIZE : int + Size of the sample batch in samples. + GPS : object + GPS band information + ''' + NAME = "Custom configuration for fast tests" + SAMPLE_RATE_HZ = freq_profile_peregrine['sampling_freq'] + SAMPLE_BATCH_SIZE = 100000 + + class GPS(object): + + class L1(object): + INTERMEDIATE_FREQUENCY_HZ = freq_profile_peregrine['GPS_L1_IF'] + INDEX = 0 + + class L2(object): + INTERMEDIATE_FREQUENCY_HZ = freq_profile_peregrine['GPS_L2_IF'] + INDEX = 1 diff --git a/peregrine/iqgen/iqgen_main.py b/peregrine/iqgen/iqgen_main.py new file mode 100644 index 0000000..009a57b --- /dev/null +++ b/peregrine/iqgen/iqgen_main.py @@ -0,0 +1,695 @@ +#!/bin/python +# Copyright (C) 2016 Swift Navigation Inc. +# Contact: Valeri Atamaniouk +# +# This source is subject to the license found in the file 'LICENSE' which must +# be be distributed together with this source. All other rights reserved. +# +# THIS CODE AND INFORMATION IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, +# EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. + +""" +The :mod:`peregrine.iqgen.iqgen_main` module contains classes and functions +related to parameter processing. + +""" +import time +import argparse +import scipy.constants +import numpy +import json +import logging +try: + import progressbar + hasProgressBar = True +except: + hasProgressBar = False + +from peregrine.iqgen.bits.satellite_gps import GPSSatellite + +# Doppler objects +from peregrine.iqgen.bits.doppler_poly import zeroDoppler +from peregrine.iqgen.bits.doppler_poly import constDoppler +from peregrine.iqgen.bits.doppler_poly import linearDoppler +from peregrine.iqgen.bits.doppler_sine import sineDoppler + +# Amplitude objects +from peregrine.iqgen.bits.amplitude_poly import AmplitudePoly +from peregrine.iqgen.bits.amplitude_sine import AmplitudeSine + +# TCXO objects +from peregrine.iqgen.bits.tcxo_poly import TCXOPoly +from peregrine.iqgen.bits.tcxo_sine import TCXOSine + +# from signals import GPS, GPS_L2C_Signal, GPS_L1CA_Signal +import peregrine.iqgen.bits.signals as signals + +from peregrine.iqgen.if_iface import LowRateConfig +from peregrine.iqgen.if_iface import NormalRateConfig +from peregrine.iqgen.if_iface import HighRateConfig +from peregrine.iqgen.if_iface import CustomRateConfig + +# Message data +from peregrine.iqgen.bits.message_const import Message as ConstMessage +from peregrine.iqgen.bits.message_zeroone import Message as ZeroOneMessage +from peregrine.iqgen.bits.message_block import Message as BlockMessage +from peregrine.iqgen.bits.message_cnav import Message as CNavMessage +from peregrine.iqgen.bits.message_lnav import Message as LNavMessage + +# PRN code generators +from peregrine.iqgen.bits.prn_gps_l1ca import PrnCode as GPS_L1CA_Code +from peregrine.iqgen.bits.prn_gps_l2c import PrnCode as GPS_L2C_Code + +# Bit stream encoders +from peregrine.iqgen.bits.encoder_gps import GPSL1BitEncoder +from peregrine.iqgen.bits.encoder_gps import GPSL2BitEncoder +from peregrine.iqgen.bits.encoder_gps import GPSL1L2BitEncoder +from peregrine.iqgen.bits.encoder_gps import GPSL1TwoBitsEncoder +from peregrine.iqgen.bits.encoder_gps import GPSL2TwoBitsEncoder +from peregrine.iqgen.bits.encoder_gps import GPSL1L2TwoBitsEncoder + +from peregrine.iqgen.generate import generateSamples + +from peregrine.iqgen.bits.satellite_factory import factoryObject as satelliteFO +from peregrine.iqgen.bits.tcxo_factory import factoryObject as tcxoFO + +logger = logging.getLogger(__name__) + + +def computeTimeDelay(doppler, symbol_index, chip_index, signal, code): + ''' + Helper function to compute signal delay to match given symbol and chip + indexes. + + Parameters + ---------- + doppler : object + Doppler object + symbol_index : long + Index of the symbol or pseudosymbol + chip_index : long + Chip index + signal : object + Signal object + code : object + Code object + + Returns + ------- + float + User's time in seconds when the user starts receiving the given symbol + and code. + ''' + if symbol_index == 0 and chip_index == 0: + return 0. + + symbolDelay_s = (1. / signal.SYMBOL_RATE_HZ) * symbol_index + chipDelay_s = (1. / signal.CODE_CHIP_RATE_HZ) * chip_index + distance_m = doppler.computeDistanceM(symbolDelay_s + chipDelay_s) + return distance_m / scipy.constants.c + + +def prepareArgsParser(): + ''' + Constructs command line argument parser. + + Returns + ------- + object + Command line argument parser object. + ''' + class AddSv(argparse.Action): + + def __init__(self, option_strings, dest, nargs=None, **kwargs): + super(AddSv, self).__init__(option_strings, dest, **kwargs) + + def __call__(self, parser, namespace, values, option_string=None): + # Initialize SV list if not yet done + if namespace.gps_sv is None: + namespace.gps_sv = [] + + # Add SV to the tail of the list. + sv = GPSSatellite(int(values)) + namespace.gps_sv.append(sv) + + # Reset all configuration parameters + namespace.l2cl_code_type = '01' + namespace.ignore_code_doppler = False + + # Doppler + namespace.doppler_type = "zero" + namespace.doppler_value = 0. + namespace.doppler_speed = 0. + namespace.distance = 0. + namespace.tec = 50. + namespace.doppler_amplitude = 0. + namespace.doppler_period = 1. + + # Source message data + namespace.message_type = "zero" + namespace.message_file = None + + # Amplitude parameters + namespace.amplitude_type = "poly" + namespace.amplitude_a0 = None + namespace.amplitude_a1 = None + namespace.amplitude_a2 = None + namespace.amplitude_a3 = None + namespace.amplitude_period = None + + class UpdateSv(argparse.Action): + + def __init__(self, option_strings, dest, nargs=None, **kwargs): + super(UpdateSv, self).__init__(option_strings, dest, **kwargs) + + def __call__(self, parser, namespace, values, option_string=None): + sv_list = getattr(namespace, "gps_sv") + if sv_list is None: + raise ValueError("No SV specified") + setattr(namespace, self.dest, values) + # super(UpdateSv, self).__call__(parser, namespace, values, option_string) + self.doUpdate(sv_list[len(sv_list) - 1], parser, namespace, values, + option_string) + + def doUpdate(self, sv, parser, namespace, values, option_string): + pass + + class UpdateBands(UpdateSv): + + def __init__(self, option_strings, dest, nargs=None, **kwargs): + super(UpdateBands, self).__init__(option_strings, dest, **kwargs) + + def doUpdate(self, sv, parser, namespace, values, option_string): + l1caEnabled = False + l2cEnabled = False + if namespace.bands == "l1ca": + l1caEnabled = True + elif namespace.bands == "l2c": + l2cEnabled = True + elif namespace.bands == "l1ca+l2c": + l1caEnabled = True + l2cEnabled = True + else: + raise ValueError() + sv.setL2CLCodeType(namespace.l2cl_code_type) + sv.setL1CAEnabled(l1caEnabled) + sv.setL2CEnabled(l2cEnabled) + + class UpdateDopplerType(UpdateSv): + + def __init__(self, option_strings, dest, nargs=None, **kwargs): + super(UpdateDopplerType, self).__init__(option_strings, dest, **kwargs) + + def doUpdate(self, sv, parser, namespace, values, option_string): + if sv.l1caEnabled: + frequency_hz = signals.GPS.L1CA.CENTER_FREQUENCY_HZ + elif sv.l2cEnabled: + frequency_hz = signals.GPS.L2C.CENTER_FREQUENCY_HZ + else: + raise ValueError("Signal band must be specified before doppler") + + if namespace.doppler_type == "zero": + doppler = zeroDoppler(namespace.distance, namespace.tec, frequency_hz) + elif namespace.doppler_type == "const": + doppler = constDoppler(namespace.distance, + namespace.tec, + frequency_hz, + namespace.doppler_value) + elif namespace.doppler_type == "linear": + doppler = linearDoppler(namespace.distance, + namespace.tec, + frequency_hz, + namespace.doppler_value, + namespace.doppler_speed) + elif namespace.doppler_type == "sine": + doppler = sineDoppler(namespace.distance, + namespace.tec, + frequency_hz, + namespace.doppler_value, + namespace.doppler_amplitude, + namespace.doppler_period) + else: + raise ValueError("Unsupported doppler type") + sv.doppler = doppler + + class DisableCodeDoppler(UpdateSv): + + def __init__(self, option_strings, dest, nargs=None, **kwargs): + super(DisableCodeDoppler, self).__init__(option_strings, dest, **kwargs) + + def doUpdate(self, sv, parser, namespace, values, option_string): + sv.getDoppler().setCodeDopplerDisabled(True) + + class UpdateAmplitudeType(UpdateSv): + + def __init__(self, option_strings, dest, nargs=None, **kwargs): + super(UpdateAmplitudeType, self).__init__(option_strings, dest, **kwargs) + + def doUpdate(self, sv, parser, namespace, values, option_string): + if namespace.amplitude_type == "poly": + coeffs = [] + hasHighOrder = False + + srcA = [namespace.amplitude_a3, namespace.amplitude_a2, + namespace.amplitude_a1, namespace.amplitude_a0] + for a in srcA: + if a is not None: + coeffs.append(a) + hasHighOrder = True + elif hasHighOrder: + coeffs.append(0.) + amplitude = AmplitudePoly(tuple(coeffs)) + elif namespace.amplitude_type == "sine": + initial = 1. + ampl = 0.5 + period_s = 1. + if namespace.amplitude_a0 is not None: + initial = namespace.amplitude_a0 + if namespace.amplitude_a1 is not None: + ampl = namespace.amplitude_a1 + if namespace.amplitude_period is not None: + period_s = namespace.amplitude_period + + amplitude = AmplitudeSine(initial, ampl, period_s) + else: + raise ValueError("Unsupported amplitude type") + sv.setAmplitude(amplitude) + + class UpdateTcxoType(UpdateSv): + + def __init__(self, option_strings, dest, nargs=None, **kwargs): + super(UpdateTcxoType, self).__init__(option_strings, dest, **kwargs) + + def doUpdate(self, sv, parser, namespace, values, option_string): + if namespace.tcxo_type == "poly": + coeffs = [] + hasHighOrder = False + + srcA = [namespace.tcxo_a3, namespace.tcxo_a2, + namespace.tcxo_a1, namespace.tcxo_a0] + for a in srcA: + if a is not None: + coeffs.append(a) + hasHighOrder = True + elif hasHighOrder: + coeffs.append(0.) + tcxo = TCXOPoly(coeffs) + elif namespace.tcxo_type == "sine": + initial = 0. + ampl = 0.5 + period_s = 1. + if namespace.tcxo_a0 is not None: + ampl = namespace.tcxo_a0 + if namespace.amplitude_a1 is not None: + ampl = namespace._a1 + if namespace.tcxo_period is not None: + period_s = namespace.tcxo_period + + tcxo = TCXOSine(initial, ampl, period_s) + else: + raise ValueError("Unsupported amplitude type") + namespace.tcxo = tcxo + + class UpdateMessageType(UpdateSv): + + def __init__(self, option_strings, dest, nargs=None, **kwargs): + super(UpdateMessageType, self).__init__(option_strings, dest, **kwargs) + + def doUpdate(self, sv, parser, namespace, values, option_string): + if namespace.message_type == "zero": + messageL1 = ConstMessage(1) + messageL2 = messageL1 + elif namespace.message_type == "one": + messageL1 = ConstMessage(-1) + messageL2 = messageL1 + elif namespace.message_type == "zero+one": + messageL1 = ZeroOneMessage() + messageL2 = messageL1 + elif namespace.message_type == "crc": + messageL1 = LNavMessage(sv.prn) + messageL2 = CNavMessage(sv.prn) + else: + raise ValueError("Unsupported message type") + sv.setL1CAMessage(messageL1) + sv.setL2CMessage(messageL2) + + class UpdateMessageFile(UpdateSv): + + def __init__(self, option_strings, dest, nargs=None, **kwargs): + super(UpdateMessageFile, self).__init__(option_strings, dest, **kwargs) + + def doUpdate(self, sv, parser, namespace, values, option_string): + data = numpy.fromfile(namespace.message_file, dtype=numpy.uint8) + namespace.message_file.close() + data = numpy.unpackbits(data) + data = numpy.asarray(data, dtype=numpy.int8) + data <<= 1 + data -= 1 + numpy.negative(data, out=data) + message = BlockMessage(data) + + sv.setL1CAMessage(message) + sv.setL2CMessage(message) + + class SaveConfigAction(argparse.Action): + + def __init__(self, option_strings, dest, nargs=None, **kwargs): + super(SaveConfigAction, self).__init__(option_strings, dest, **kwargs) + + def __call__(self, parser, namespace, values, option_string=None): + + gps_sv = namespace.gps_sv + + encoded_gps_sv = [satelliteFO.toMapForm(sv) for sv in gps_sv] + + data = {'type': 'Namespace', + 'gps_sv': encoded_gps_sv, + 'profile': namespace.profile, + 'encoder': namespace.encoder, + 'chip_delay': namespace.chip_delay, + 'symbol_delay': namespace.symbol_delay, + 'generate': namespace.generate, + 'snr': namespace.snr, + 'filter_type': namespace.filter_type, + 'tcxo': tcxoFO.toMapForm(namespace.tcxo) + } + json.dump(data, values, indent=2) + values.close() + namespace.no_run = True + + class LoadConfigAction(argparse.Action): + + def __init__(self, option_strings, dest, nargs=None, **kwargs): + super(LoadConfigAction, self).__init__(option_strings, dest, **kwargs) + + def __call__(self, parser, namespace, values, option_string=None): + loaded = json.load(values) + namespace.profile = loaded['profile'] + namespace.encoder = loaded['encoder'] + namespace.chip_delay = loaded['chip_delay'] + namespace.symbol_delay = loaded['symbol_delay'] + namespace.generate = loaded['generate'] + namespace.snr = loaded['snr'] + namespace.filter_type = loaded['filter_type'] + namespace.tcxo = tcxoFO.fromMapForm(loaded['tcxo']) + namespace.gps_sv = [ + satelliteFO.fromMapForm(sv) for sv in loaded['gps_sv']] + values.close() + + parser = argparse.ArgumentParser( + description="Signal generator", usage='%(prog)s [options]') + parser.add_argument('--gps-sv', + default=[], + help='Enable GPS satellite', + action=AddSv) + parser.add_argument('--bands', + default="l1ca", + choices=["l1ca", "l2c", "l1ca+l2c"], + help="Signal bands to enable", + action=UpdateBands) + parser.add_argument('--l2cl-code-type', + default='01', + choices=['01', '1', '0'], + help="GPS L2 CL code type", + action=UpdateBands) + parser.add_argument('--doppler-type', + default="zero", + choices=["zero", "const", "linear", "sine"], + help="Configure doppler type", + action=UpdateDopplerType) + parser.add_argument('--doppler-value', + type=float, + help="Doppler shift in hertz (initial)", + action=UpdateDopplerType) + parser.add_argument('--doppler-speed', + type=float, + help="Doppler shift change in hertz/second", + action=UpdateDopplerType) + parser.add_argument('--distance', + type=float, + help="Distance in meters for signal delay (initial)", + action=UpdateDopplerType) + parser.add_argument('--tec', + type=float, + help="Ionosphere TEC for signal delay" + " (electrons per meter^2)", + action=UpdateDopplerType) + parser.add_argument('--doppler-amplitude', + type=float, + help="Doppler change amplitude (hertz)", + action=UpdateDopplerType) + parser.add_argument('--doppler-period', + type=float, + help="Doppler change period (seconds)", + action=UpdateDopplerType) + parser.add_argument('--ignore-code-doppler', + help="Disable doppler for code and data processing", + action=DisableCodeDoppler) + parser.add_argument('--amplitude-type', + default="poly", + choices=["poly", "sine"], + help="Configure amplitude type: polynomial or sine.", + action=UpdateAmplitudeType) + parser.add_argument('--amplitude-a0', + type=float, + help="Amplitude coefficient (a0 for polynomial;" + " offset for sine)", + action=UpdateAmplitudeType) + parser.add_argument('--amplitude-a1', + type=float, + help="Amplitude coefficient (a1 for polynomial," + " amplitude for size)", + action=UpdateAmplitudeType) + parser.add_argument('--amplitude-a2', + type=float, + help="Amplitude coefficient (a2 for polynomial)", + action=UpdateAmplitudeType) + parser.add_argument('--amplitude-a3', + type=float, + help="Amplitude coefficient (a3 for polynomial)", + action=UpdateAmplitudeType) + parser.add_argument('--amplitude-period', + type=float, + help="Amplitude period in seconds for sine", + action=UpdateAmplitudeType) + parser.add_argument('--message-type', default="zero", + choices=["zero", "one", "zero+one", "crc"], + help="Message type", + action=UpdateMessageType) + parser.add_argument('--message-file', + type=argparse.FileType('rb'), + help="Source file for message contents.", + action=UpdateMessageFile) + parser.add_argument('--symbol_delay', + type=int, + help="Initial symbol index") + parser.add_argument('--chip_delay', + type=int, + help="Initial chip index") + parser.add_argument('--filter-type', + default='none', + choices=['none', 'lowpass', 'bandpass'], + help="Enable filter") + parser.add_argument('--snr', + type=float, + help="SNR for noise generation") + parser.add_argument('--tcxo-type', + choices=["poly", "sine"], + help="TCXO drift type", + action=UpdateTcxoType) + parser.add_argument('--tcxo-a0', + type=float, + help="TCXO a0 coefficient for polynomial TCXO drift" + " or initial shift for sine TCXO drift", + action=UpdateTcxoType) + parser.add_argument('--tcxo-a1', + type=float, + help="TCXO a1 coefficient for polynomial TCXO drift" + " or amplitude for sine TCXO drift", + action=UpdateTcxoType) + parser.add_argument('--tcxo-a2', + type=float, + help="TCXO a2 coefficient for polynomial TCXO drift", + action=UpdateTcxoType) + parser.add_argument('--tcxo-a3', + type=float, + help="TCXO a3 coefficient for polynomial TCXO drift", + action=UpdateTcxoType) + parser.add_argument('--tcxo-period', + type=float, + help="TCXO period in seconds for sine TCXO drift", + action=UpdateTcxoType) + parser.add_argument('--debug', + type=argparse.FileType('wb'), + help="Debug output file") + parser.add_argument('--generate', + type=float, + default=3., + help="Amount of data to generate, in seconds") + parser.add_argument('--encoder', + default="2bits", + choices=["1bit", "2bits"], + help="Output data format") + parser.add_argument('--output', + type=argparse.FileType('wb'), + help="Output file name") + parser.add_argument('--profile', + default="normal_rate", + choices=["low_rate", "normal_rate", "high_rate", + "custom_rate"], + help="Output profile configuration") + parser.add_argument('-j', '--jobs', + type=int, + default=0, + help="Use parallel threads") + + parser.add_argument('--save-config', + type=argparse.FileType('wt'), + help="Store configuration into file (implies --no-run)", + action=SaveConfigAction) + + parser.add_argument('--load-config', + type=argparse.FileType('rt'), + help="Restore configuration from file", + action=LoadConfigAction) + + parser.add_argument('--no-run', + action="store_true", + default=False, + help="Do not generate output.") + + parser.set_defaults(tcxo=TCXOPoly(())) + + return parser + + +def main(): + from peregrine.log import default_logging_config + default_logging_config() + + parser = prepareArgsParser() + args = parser.parse_args() + + if args.no_run: + return 0 + + if args.output is None: + parser.print_help() + return 0 + + if args.profile == "low_rate": + outputConfig = LowRateConfig + elif args.profile == "normal_rate": + outputConfig = NormalRateConfig + elif args.profile == "high_rate": + outputConfig = HighRateConfig + elif args.profile == "custom_rate": + outputConfig = CustomRateConfig + else: + raise ValueError() + + print "Output configuration:" + print " Description: ", outputConfig.NAME + print " Sampling rate: ", outputConfig.SAMPLE_RATE_HZ + print " Batch size: ", outputConfig.SAMPLE_BATCH_SIZE + print " GPS L1 IF: ", outputConfig.GPS.L1.INTERMEDIATE_FREQUENCY_HZ + print " GPS L2 IF: ", outputConfig.GPS.L2.INTERMEDIATE_FREQUENCY_HZ + print "Other parameters:" + print " TCXO: ", args.tcxo + print " SNR: ", args.snr + print " tSatellites: ", args.gps_sv + + # Check which signals are enabled on each of satellite to select proper + # output encoder + enabledGPSL1 = False + enabledGPSL2 = False + + for sv in args.gps_sv: + enabledGPSL1 |= sv.isBandEnabled(outputConfig.GPS.L1.INDEX, outputConfig) + enabledGPSL2 |= sv.isBandEnabled(outputConfig.GPS.L2.INDEX, outputConfig) + + # Configure data encoder + if args.encoder == "1bit": + if enabledGPSL1 and enabledGPSL2: + encoder = GPSL1L2BitEncoder(outputConfig) + elif enabledGPSL2: + encoder = GPSL2BitEncoder(outputConfig) + else: + encoder = GPSL1BitEncoder(outputConfig) + elif args.encoder == "2bits": + if enabledGPSL1 and enabledGPSL2: + encoder = GPSL1L2TwoBitsEncoder(outputConfig) + elif enabledGPSL2: + encoder = GPSL2TwoBitsEncoder(outputConfig) + else: + encoder = GPSL1TwoBitsEncoder(outputConfig) + else: + raise ValueError("Encoder type is not supported") + + if enabledGPSL1: + signal = signals.GPS.L1CA + code = GPS_L1CA_Code + elif enabledGPSL2: + signal = signals.GPS.L2C + code = GPS_L2C_Code + else: + signal = signals.GPS.L1CA + code = GPS_L1CA_Code + + # Compute time delay for the needed bit/chip number + # This delay is computed for the first satellite + initial_symbol_idx = 0 # Initial symbol index + initial_chip_idx = 0 # Initial chip index + if args.chip_delay is not None: + initial_chip_idx = args.chip_delay + if args.symbol_delay is not None: + initial_chip_idx = args.symbol_delay + + time0_s = computeTimeDelay(args.gps_sv[0].doppler, + initial_symbol_idx, + initial_chip_idx, + signal, + code) + logger.debug("Computed symbol/chip delay={} seconds".format(time0_s)) + + startTime_s = time.time() + n_samples = long(outputConfig.SAMPLE_RATE_HZ * args.generate) + + logger.debug("Generating {} samples for {} seconds". + format(n_samples, args.generate)) + + if hasProgressBar: + widgets = ['Generating ', + progressbar.Counter(), ' ', + progressbar.Percentage(), ' ', + progressbar.ETA(), ' ', + progressbar.Bar()] + pbar = progressbar.ProgressBar(widgets=widgets, + maxval=n_samples).start() + else: + pbar = None + + generateSamples(args.output, + args.gps_sv, + encoder, + time0_s, + n_samples, + outputConfig, + tcxo=args.tcxo, + SNR=args.snr, + filterType=args.filter_type, + logFile=args.debug, + threadCount=args.jobs, + pbar=pbar) + args.output.close() + # if pbar: + # pbar.finish() + + duration_s = time.time() - startTime_s + ratio = n_samples / duration_s + logger.debug("Total time = {} sec. Ratio={} samples per second". + format(duration_s, ratio)) + +if __name__ == '__main__': + main() diff --git a/peregrine/parallel_processing.py b/peregrine/parallel_processing.py index 44238bb..be8b79b 100644 --- a/peregrine/parallel_processing.py +++ b/peregrine/parallel_processing.py @@ -4,70 +4,83 @@ import progressbar as pb import multiprocessing as mp import time +import sys +import traceback + def spawn(f): - def worker(q_in, q_out, q_progress): - while True: - i,x = q_in.get() - if i is None: - break - try: - if q_progress: - q_out.put((i,f(x, q_progress=q_progress))) - else: - q_out.put((i,f(x))) - except Exception as err: - print "Subprocess raised exception:" - print err - q_out.put(None) - return worker + def worker(q_in, q_out, q_progress): + while True: + i, x = q_in.get() + if i is None: + break + try: + if q_progress: + res = f(x, q_progress=q_progress) + q_out.put((i, res)) + else: + res = f(x) + q_out.put((i, res)) + except: + print "Subprocess raised exception:" + exType, exValue, exTraceback = sys.exc_info() + traceback.print_exception( + exType, exValue, exTraceback, file=sys.stdout) + q_out.put(None) + return worker -def parmap(f, X, nprocs = mp.cpu_count(), show_progress=True, func_progress=False): - q_in = mp.Queue() - q_out = mp.Queue() - if func_progress: - q_progress = mp.Queue(100) - else: - q_progress = None - proc = [mp.Process(target=spawn(f),args=(q_in, q_out, q_progress)) for _ in range(nprocs)] +def parmap(f, X, nprocs=mp.cpu_count(), show_progress=True, func_progress=False): + q_in = mp.Queue() + q_out = mp.Queue() + if func_progress: + q_progress = mp.Queue(100) + else: + q_progress = None - for p in proc: - p.daemon = True - p.start() + if nprocs > mp.cpu_count(): + nprocs = mp.cpu_count() - if show_progress: - pbar = pb.ProgressBar(widgets=[pb.Percentage(), ' ', pb.ETA()], maxval=len(X)).start() + proc = [mp.Process(target=spawn(f), args=(q_in, q_out, q_progress)) + for _ in range(nprocs)] - [q_in.put((i, x)) for i, x in enumerate(X)] + for p in proc: + p.daemon = True + p.start() - [q_in.put((None,None)) for _ in range(nprocs)] + if show_progress: + pbar = pb.ProgressBar( + widgets=[pb.Percentage(), ' ', pb.ETA()], maxval=len(X)).start() - n_done = 0 - progress = 0 - res = [] - t0 = time.time() - while n_done < len(X): - if func_progress: - time.sleep(0.02) - else: - res.append(q_out.get()) - n_done += 1 - while not q_out.empty(): - res.append(q_out.get()) - n_done += 1 - if q_progress: - while not q_progress.empty(): - progress_increment = q_progress.get_nowait() - progress += progress_increment - else: - progress = n_done - if show_progress and progress <= len(X): - pbar.update(progress) + [q_in.put((i, x)) for i, x in enumerate(X)] + + [q_in.put((None, None)) for _ in range(nprocs)] + + n_done = 0 + progress = 0 + res = [] + t0 = time.time() + while n_done < len(X): + if func_progress: + time.sleep(0.02) + else: + res.append(q_out.get()) + n_done += 1 + while not q_out.empty(): + res.append(q_out.get()) + n_done += 1 + if q_progress: + while not q_progress.empty(): + progress_increment = q_progress.get_nowait() + progress += progress_increment + else: + progress = n_done + if show_progress and progress <= len(X): + pbar.update(progress) - if show_progress: - pbar.finish() + if show_progress: + pbar.finish() - [p.join() for p in proc] + [p.join() for p in proc] - return [x for i,x in sorted(res)] + return [x for i, x in sorted(res)] diff --git a/peregrine/run.py b/peregrine/run.py index 7e7eb86..a0c6500 100755 --- a/peregrine/run.py +++ b/peregrine/run.py @@ -13,27 +13,44 @@ import argparse import cPickle import logging -from operator import attrgetter import numpy as np +from operator import attrgetter from peregrine.samples import load_samples -from peregrine.acquisition import Acquisition, load_acq_results, save_acq_results +from peregrine.acquisition import Acquisition, load_acq_results,\ + save_acq_results from peregrine.navigation import navigation -from peregrine.tracking import track +import peregrine.tracking as tracking from peregrine.log import default_logging_config -import defaults +from peregrine import defaults +from peregrine.initSettings import initSettings +from peregrine.gps_constants import L1CA, L2C, GLO_L1 + + +def unpickle_iter(filenames): + try: + f = [open(filename, "r") for filename in filenames] + + while True: + yield [cPickle.load(fh) for fh in f] + + except EOFError: + raise StopIteration + + finally: + for fh in f: + fh.close() -from initSettings import initSettings def main(): default_logging_config() - # Initialize constants, settings - settings = initSettings() - parser = argparse.ArgumentParser() parser.add_argument("file", help="the sample data file to process") + parser.add_argument("--skip-glonass", + help="skip glonass", + action="store_true") parser.add_argument("-a", "--skip-acquisition", help="use previously saved acquisition results", action="store_true") @@ -43,15 +60,101 @@ def main(): parser.add_argument("-n", "--skip-navigation", help="use previously saved navigation results", action="store_true") + parser.add_argument("--ms-to-process", + help="the number of milliseconds to process." + "(-1: use all available data", + default="-1") + parser.add_argument("--profile", + help="L1C/A & L2C IF + sampling frequency profile" + "('peregrine'/'custom_rate', 'low_rate', " + "'normal_rate' (piksi_v3), 'high_rate')", + default='peregrine') parser.add_argument("-f", "--file-format", default=defaults.file_format, help="the format of the sample data file " - "(e.g. 'piksi', 'int8', '1bit', '1bitrev')") + "('piksi', 'int8', '1bit', '1bitrev', " + "'1bit_x2', '2bits', '2bits_x2', '2bits_x4')") + parser.add_argument('--l1ca-profile', + help='L1 C/A stage profile', + choices=defaults.l1ca_stage_profiles.keys()) + if sys.stdout.isatty(): + progress_bar_default = 'stdout' + elif sys.stderr.isatty(): + progress_bar_default = 'stderr' + else: + progress_bar_default = 'none' + parser.add_argument("--progress-bar", + metavar='FLAG', + choices=['stdout', 'stderr', 'none'], + default=progress_bar_default, + help="Show progress bar. Default is '%s'" % + progress_bar_default) + + fpgaSim = parser.add_argument_group('FPGA simulation', + 'FPGA delay control simulation') + + fpgaSim.add_argument("--pipelining", + type=float, + nargs='?', + help="Use FPGA pipelining simulation. Supply optional " + " coefficient (%f)" % defaults.pipelining_k, + const=defaults.pipelining_k, + default=None) + + fpgaSim.add_argument("--short-long-cycles", + type=float, + nargs='?', + help="Use FPGA short-long cycle simulation. Supply" + " optional pipelining coefficient (0.)", + const=0., + default=None) args = parser.parse_args() + + if args.file is None: + parser.print_help() + return + + if args.profile == 'peregrine' or args.profile == 'custom_rate': + freq_profile = defaults.freq_profile_peregrine + elif args.profile == 'low_rate': + freq_profile = defaults.freq_profile_low_rate + elif args.profile == 'normal_rate': + freq_profile = defaults.freq_profile_normal_rate + elif args.profile == 'high_rate': + freq_profile = defaults.freq_profile_high_rate + else: + raise NotImplementedError() + + if args.l1ca_profile: + profile = defaults.l1ca_stage_profiles[args.l1ca_profile] + stage2_coherent_ms = profile[1]['coherent_ms'] + stage2_params = profile[1]['loop_filter_params'] + else: + stage2_coherent_ms = None + stage2_params = None + + if args.pipelining is not None: + tracker_options = {'mode': 'pipelining', 'k': args.pipelining} + else: + tracker_options = None + + settings = initSettings(freq_profile) settings.fileName = args.file + ms_to_process = int(args.ms_to_process) + samplesPerCode = int(round(settings.samplingFreq / (settings.codeFreqBasis / settings.codeLength))) + gloSamplesPerCode = int(round(settings.samplingFreq / + (settings.gloCodeFreqBasis / + settings.gloCodeLength))) + + samples = {L1CA: {'IF': freq_profile['GPS_L1_IF']}, + L2C: {'IF': freq_profile['GPS_L2_IF']}, + GLO_L1: {'IF': freq_profile['GLO_L1_IF']}, + 'samples_total': -1, + 'sample_index': settings.skipNumberOfBytes} + # Do acquisition acq_results_file = args.file + ".acq_results" if args.skip_acquisition: @@ -64,11 +167,28 @@ def main(): sys.exit(1) else: # Get 11ms of acquisition samples for fine frequency estimation - acq_samples = load_samples(args.file, 11*samplesPerCode, - settings.skipNumberOfBytes, - file_format=args.file_format) - acq = Acquisition(acq_samples) - acq_results = acq.acquisition() + load_samples(samples=samples, + num_samples=11 * max(samplesPerCode, gloSamplesPerCode), + filename=args.file, + file_format=args.file_format) + + acq = Acquisition(samples[L1CA]['samples'], + freq_profile['sampling_freq'], + freq_profile['GPS_L1_IF'], + defaults.code_period * freq_profile['sampling_freq']) + acq_results = acq.acquisition(progress_bar_output=args.progress_bar) + + if not args.skip_glonass: + acq = Acquisition(samples[GLO_L1]['samples'], + freq_profile['sampling_freq'], + freq_profile['GLO_L1_IF'], + defaults.glo_code_period * + freq_profile['sampling_freq'], + code_length=settings.gloCodeLength) + acq_results = acq.acquisition(bandcode=GLO_L1, + progress_bar_output=args.progress_bar) + + print "Acquisition is over!" try: save_acq_results(acq_results_file, acq_results) @@ -98,36 +218,60 @@ def main(): track_results_file) sys.exit(1) else: - signal = load_samples(args.file, - int(settings.samplingFreq*1e-3*(settings.msToProcess+22)), - settings.skipNumberOfBytes, - file_format=args.file_format) - track_results = track(signal, acq_results, settings.msToProcess) - try: - with open(track_results_file, 'wb') as f: - cPickle.dump(track_results, f, protocol=cPickle.HIGHEST_PROTOCOL) - logging.debug("Saving tracking results as '%s'" % track_results_file) - except IOError: - logging.error("Couldn't save tracking results file '%s'.", - track_results_file) + load_samples(samples=samples, + filename=args.file, + file_format=args.file_format) + + if ms_to_process < 0: + ms_to_process = int( + 1e3 * samples['samples_total'] / freq_profile['sampling_freq']) + + tracker = tracking.Tracker(samples=samples, + channels=acq_results, + ms_to_track=ms_to_process, + sampling_freq=freq_profile[ + 'sampling_freq'], # [Hz] + stage2_coherent_ms=stage2_coherent_ms, + stage2_loop_filter_params=stage2_params, + tracker_options=tracker_options, + output_file=args.file, + progress_bar_output=args.progress_bar) + tracker.start() + condition = True + while condition: + sample_index = tracker.run_channels(samples) + if sample_index == samples['sample_index']: + condition = False + else: + samples['sample_index'] = sample_index + load_samples(samples=samples, + filename=args.file, + file_format=args.file_format) + fn_results = tracker.stop() + + logging.debug("Saving tracking results as '%s'" % fn_results) # Do navigation nav_results_file = args.file + ".nav_results" if not args.skip_navigation: - nav_solns = navigation(track_results, settings) - nav_results = [] - for s, t in nav_solns: - nav_results += [(t, s.pos_llh, s.vel_ned)] - if len(nav_results): - print "First nav solution: t=%s lat=%.5f lon=%.5f h=%.1f vel_ned=(%.2f, %.2f, %.2f)" % ( - nav_results[0][0], - np.degrees(nav_results[0][1][0]), np.degrees(nav_results[0][1][1]), nav_results[0][1][2], - nav_results[0][2][0], nav_results[0][2][1], nav_results[0][2][2]) - with open(nav_results_file, 'wb') as f: - cPickle.dump(nav_results, f, protocol=cPickle.HIGHEST_PROTOCOL) - print "and %d more are cPickled in '%s'." % (len(nav_results)-1, nav_results_file) - else: - print "No navigation results." + track_results_generator = lambda: unpickle_iter(fn_results) + for track_results in track_results_generator(): + nav_solns = navigation(track_results_generator, settings) + nav_results = [] + for s, t in nav_solns: + nav_results += [(t, s.pos_llh, s.vel_ned)] + if len(nav_results): + print "First nav solution: t=%s lat=%.5f lon=%.5f h=%.1f vel_ned=(%.2f, %.2f, %.2f)" % ( + nav_results[0][0], + np.degrees(nav_results[0][1][0]), np.degrees( + nav_results[0][1][1]), nav_results[0][1][2], + nav_results[0][2][0], nav_results[0][2][1], nav_results[0][2][2]) + with open(nav_results_file, 'wb') as f: + cPickle.dump(nav_results, f, protocol=cPickle.HIGHEST_PROTOCOL) + print "and %d more are cPickled in '%s'." % (len(nav_results) - 1, + nav_results_file) + else: + print "No navigation results." if __name__ == '__main__': main() diff --git a/peregrine/samples.py b/peregrine/samples.py index cc5fefc..558edb3 100644 --- a/peregrine/samples.py +++ b/peregrine/samples.py @@ -9,11 +9,140 @@ """Functions for handling sample data and sample data files.""" +import os import numpy as np +import math +import defaults +from peregrine.gps_constants import L1CA, L2C, GLO_L1 __all__ = ['load_samples', 'save_samples'] -def load_samples(filename, num_samples=-1, num_skip=0, file_format='piksi'): +def __load_samples_n_bits(filename, num_samples, num_skip, n_rx, n_bits, + value_lookup, channel_lookup = None): + ''' + Helper method to load two-bit samples from a file. + + Parameters + ---------- + filename : string + Filename of sample data file. + num_samples : int + Number of samples to read, ``-1`` means the whole file. + num_skip : int + Number of samples to discard from the beginning of the file. + n_rx : int + Number of interleaved streams in the source file + n_bits : int + Number of bits per sample + channel_lookup : array-like + Array to map channels + value_lookup : array-like + Array to map values + + Returns + ------- + out : :class:`numpy.ndarray`, shape(`n_rx`, `num_samples`,) + The sample data as a two-dimensional numpy array. The first dimension + separates codes (bands). The second dimention contains samples indexed + with the `value_lookup` table. + ''' + if not channel_lookup: + channel_lookup = range(n_rx) + + sample_block_size = n_bits * n_rx + byte_offset = int(math.floor(num_skip / (8 / sample_block_size))) + sample_offset = num_skip % (8 / sample_block_size) + s_file = np.memmap(filename, offset=byte_offset, dtype=np.uint8, mode='r') + + if num_samples > 0: + # Number of samples is defined: trim the source from end + s_file = s_file[:(num_samples * sample_block_size + 7) / 8] + + num_samples = len(s_file) * 8 / sample_block_size + + # Compute total data block size to ignore bits in the tail. + rounded_len = num_samples * sample_block_size + + bits = np.unpackbits(s_file) + samples = np.empty((n_rx, num_samples - sample_offset), + dtype=value_lookup.dtype) + + for rx in range(n_rx): + # Construct multi-bit sample values + tmp = bits[rx * n_bits:rounded_len:sample_block_size] + for bit in range(1, n_bits): + tmp <<= 1 + tmp += bits[rx * n_bits + bit:rounded_len:sample_block_size] + # Generate sample values using value_lookup table + chan = value_lookup[tmp] + chan = chan[sample_offset:] + samples[channel_lookup[rx]][:] = chan + return samples + +def __load_samples_one_bit(filename, num_samples, num_skip, n_rx, + channel_lookup = None): + ''' + Helper method to load single-bit samples from a file. + + Parameters + ---------- + filename : string + Filename of sample data file. + num_samples : int + Number of samples to read, ``-1`` means the whole file. + num_skip : int + Number of samples to discard from the beginning of the file. + n_rx : int + Number of interleaved streams in the source file + channel_lookup : array-like + Array to map channels + + Returns + ------- + out : :class:`numpy.ndarray`, shape(`n_rx`, `num_samples`,) + The sample data as a two-dimensional numpy array. The first dimension + separates codes (bands). The second dimention contains samples with one + of the values: -1, 1 + ''' + value_lookup = np.asarray((1, -1), dtype=np.int8) + return __load_samples_n_bits(filename, num_samples, num_skip, n_rx, 1, + value_lookup, channel_lookup) + +def __load_samples_two_bits(filename, num_samples, num_skip, n_rx, + channel_lookup = None): + ''' + Helper method to load two-bit samples from a file. + + Parameters + ---------- + filename : string + Filename of sample data file. + num_samples : int + Number of samples to read, ``-1`` means the whole file. + num_skip : int + Number of samples to discard from the beginning of the file. + n_rx : int + Number of interleaved streams in the source file + channel_lookup : array-like + Array to map channels + + Returns + ------- + out : :class:`numpy.ndarray`, shape(`n_rx`, `num_samples`,) + The sample data as a two-dimensional numpy array. The first dimension + separates codes (bands). The second dimention contains samples with one + of the values: -3, -1, 1, 3 + ''' + # Interleaved two bit samples from two receivers. First bit is a sign of the + # sample, and the second bit is the amplitude value: 1 or 3. + value_lookup = np.asarray((-1, -3, 1, 3), dtype=np.int8) + return __load_samples_n_bits(filename, num_samples, num_skip, n_rx, 2, + value_lookup, channel_lookup) + +def _load_samples(filename, + num_samples = defaults.processing_block_size, + num_skip = 0, + file_format = 'piksi'): """ Load sample data from a file. @@ -42,8 +171,9 @@ def load_samples(filename, num_samples=-1, num_skip=0, file_format='piksi'): Returns ------- - out : :class:`numpy.ndarray`, shape(`num_samples`,) - The sample data as a numpy array. + out : :class:`numpy.ndarray`, shape(bands, `num_samples`,) + The sample data as a two-dimensional numpy array. The first dimension + separates codes (bands). Raises ------ @@ -58,37 +188,38 @@ def load_samples(filename, num_samples=-1, num_skip=0, file_format='piksi'): if file_format == 'int8': with open(filename, 'rb') as f: f.seek(num_skip) - samples = np.fromfile(f, dtype=np.int8, count=num_samples) + samples = np.zeros((1, num_samples), dtype=np.int8) + samples[:] = np.fromfile(f, dtype=np.int8, count=num_samples) elif file_format == 'c8c8': # Interleaved complex samples from two receivers, i.e. first four bytes are # I0 Q0 I1 Q1 s_file = np.memmap(filename, offset=num_skip, dtype=np.int8, mode='r') n_rx = 2 if num_samples > 0: - s_file = s_file[:num_samples*2*n_rx] - samples = np.empty([n_rx, len(s_file)/(2 * n_rx)], dtype=np.complex64) + s_file = s_file[:num_samples * 2 * n_rx] + samples = np.empty([n_rx, len(s_file) / (2 * n_rx)], dtype=np.complex64) for rx in range(n_rx): - samples[rx] = s_file[2 * rx : : 2 * n_rx] + s_file[2*rx + 1 :: 2 * n_rx]*1j + samples[rx] = s_file[2 * rx : : 2 * n_rx] + s_file[2 * rx + 1 :: 2 * n_rx] * 1j elif file_format == 'c8c8_tayloe': # Interleaved complex samples from two receivers, i.e. first four bytes are # I0 Q0 I1 Q1. Tayloe-upconverted to become purely real with fs=4fs0, fi=fs0 s_file = np.memmap(filename, offset=num_skip, dtype=np.int8, mode='r') n_rx = 2 if num_samples > 0: - s_file = s_file[:num_samples*2*n_rx] - samples = np.empty([n_rx, 4*len(s_file)/(2 * n_rx)], dtype=np.int8) + s_file = s_file[:num_samples * 2 * n_rx] + samples = np.empty([n_rx, 4 * len(s_file) / (2 * n_rx)], dtype=np.int8) for rx in range(n_rx): - samples[rx][0::4] = s_file[2 * rx : : 2 * n_rx] + samples[rx][0::4] = s_file[2 * rx : : 2 * n_rx] samples[rx][1::4] = -s_file[2 * rx + 1 : : 2 * n_rx] samples[rx][2::4] = -s_file[2 * rx : : 2 * n_rx] - samples[rx][3::4] = s_file[2 * rx + 1 : : 2 * n_rx] + samples[rx][3::4] = s_file[2 * rx + 1 : : 2 * n_rx] elif file_format == 'piksinew': packed = np.memmap(filename, offset=num_skip, dtype=np.uint8, mode='r') if num_samples > 0: packed = packed[:num_samples] - samples = np.empty(len(packed), dtype=np.int8) - samples[:] = (packed >> 6) - 1 + samples = np.empty((1, len(packed)), dtype=np.int8) + samples[0][:] = (packed >> 6) - 1 elif file_format == 'piksi': """ @@ -105,7 +236,7 @@ def load_samples(filename, num_samples=-1, num_skip=0, file_format='piksi'): if num_samples > 0: num_skip_bytes = num_skip / 2 num_skip_samples = num_skip % 2 - num_bytes = num_samples/2 + 1 + num_bytes = num_samples / 2 + 1 else: # Read whole file num_skip_bytes = 0 @@ -122,17 +253,20 @@ def load_samples(filename, num_samples=-1, num_skip=0, file_format='piksi'): samples[::2] = (packed >> 5) samples[1::2] = (packed >> 2) & 7 # Sign-magnitude to two's complement mapping - samples = (1-2*(samples>>2)) * (2*(samples&3)+1) + samples = (1 - 2 * (samples >> 2)) * (2 * (samples & 3) + 1) samples = samples[num_skip_samples:] if num_samples > 0: samples = samples[:num_samples] + tmp = np.ndarray((1, len(samples)), dtype=np.int8) + tmp[0][:] = samples + samples = tmp elif file_format == '1bit' or file_format == '1bitrev': if num_samples > 0: num_skip_bytes = num_skip / 8 num_skip_samples = num_skip % 8 - num_bytes = num_samples/8 + 1 + num_bytes = num_samples / 8 + 1 else: # Read whole file num_skip_bytes = 0 @@ -146,17 +280,93 @@ def load_samples(filename, num_samples=-1, num_skip=0, file_format='piksi'): samples *= 2 samples -= 1 if file_format == '1bitrev': - samples = np.reshape(samples,(-1,8))[:,::-1].flatten(); + samples = np.reshape(samples, (-1, 8))[:, ::-1].flatten() samples = samples[num_skip_samples:] if num_samples > 0: samples = samples[:num_samples] + tmp = np.ndarray((1, len(samples)), dtype=np.int8) + tmp[0][:] = samples + samples = tmp + elif file_format == '1bit_x2': + # Interleaved single bit samples from two receivers: -1, +1 + samples = __load_samples_one_bit(filename, num_samples, num_skip, 2, + defaults.file_encoding_1bit_x2) + elif file_format == '2bits': + # Two bit samples from one receiver: -3, -1, +1, +3 + samples = __load_samples_two_bits(filename, num_samples, num_skip, 1) + elif file_format == '2bits_x2': + # Interleaved two bit samples from two receivers: -3, -1, +1, +3 + samples = __load_samples_two_bits(filename, num_samples, num_skip, 2, + defaults.file_encoding_2bits_x2) + elif file_format == '2bits_x4': + # Interleaved two bit samples from four receivers: -3, -1, +1, +3 + samples = __load_samples_two_bits(filename, num_samples, num_skip, 4, + defaults.file_encoding_2bits_x4) else: raise ValueError("Unknown file type '%s'" % file_format) - if len(samples.T) < num_samples: - raise EOFError("Failed to read %d samples from file '%s'" % - (num_samples, filename)) + return samples + +def __get_samples_total(filename, file_format, sample_index): + if file_format == 'int8': + samples_block_size = 8 + elif file_format == 'piksi': + """ + Piksi format is packed 3-bit sign-magnitude samples, 2 samples per byte. + + Bits: + [1..0] Flags (reserved for future use) + [3..2] Sample 2 magnitude + [4] Sample 2 sign (1 is -ve) + [6..5] Sample 1 magnitude + [7] Sample 1 sign (1 is -ve) + + """ + samples_block_size = 4 + elif file_format == '1bit' or file_format == '1bitrev': + samples_block_size = 1 + elif file_format == '1bit_x2': + # Interleaved single bit samples from two receivers: -1, +1 + samples_block_size = 2 + elif file_format == '2bits': + # Two bit samples from one receiver: -3, -1, +1, +3 + samples_block_size = 2 + elif file_format == '2bits_x2': + # Interleaved two bit samples from two receivers: -3, -1, +1, +3 + samples_block_size = 4 + elif file_format == '2bits_x4': + # Interleaved two bit samples from four receivers: -3, -1, +1, +3 + samples_block_size = 8 + else: + raise ValueError("Unknown file type '%s'" % file_format) + + file_size = os.path.getsize(filename) + samples_total = 8 * file_size / samples_block_size + + if sample_index < samples_total: + samples_total -= sample_index + + return samples_total + +def load_samples(samples, + filename, + num_samples = defaults.processing_block_size, + file_format = 'piksi'): + + if samples['samples_total'] == -1: + samples['samples_total'] = __get_samples_total(filename, + file_format, + samples['sample_index']) + signal = _load_samples(filename, + num_samples, + samples['sample_index'], + file_format) + samples[L1CA]['samples'] = signal[defaults.sample_channel_GPS_L1] + if len(signal) > 1: + samples[L2C]['samples'] = signal[defaults.sample_channel_GPS_L2] + if len(signal) > 2: + samples[GLO_L1]['samples'] = signal[defaults.sample_channel_GLO_L1] return samples diff --git a/peregrine/stream_usrp.py b/peregrine/stream_usrp.py index 9abeae3..29832c5 100755 --- a/peregrine/stream_usrp.py +++ b/peregrine/stream_usrp.py @@ -10,7 +10,6 @@ from gnuradio import blocks from gnuradio import gr from gnuradio import uhd -from gnuradio import gru import time import argparse import sys @@ -24,6 +23,7 @@ def __init__(self, filenames, dev_addrs, dual, gr.top_block.__init__(self) if mix: raise NotImplementedError("TODO: Hilbert remix mode not implemented.") + if dual: channels = [0, 1] else: @@ -58,6 +58,7 @@ def __init__(self, filenames, dev_addrs, dual, if noise or onebit or not iq: raise NotImplementedError("TODO: RX channel-interleaved mode only " "supported for noiseless 8-bit complex.") + BLOCK_N = 16*1024*1024 demux = blocks.vector_to_streams(2, len(uhd_sinks)) self.connect(blocks.file_source(2*len(uhd_sinks)*BLOCK_N, filenames[0], False), @@ -131,23 +132,49 @@ def __init__(self, filenames, dev_addrs, dual, # No external PPS/10 MHz. Just set each clock and accept some skew. t = time.time() [sink.set_time_now(uhd.time_spec(time.time())) for sink in uhd_sinks] - if len(uhd_sinks) > 1 or dual: + if len(uhd_sinks) > 1: print "Uncabled; loosely synced only. Initial skew ~ %.1f ms" % ( (time.time()-t) * 1000) t_start = uhd.time_spec(time.time() + 1.5) [sink.set_start_time(t_start) for sink in uhd_sinks] print "ready" - # setup message handler - self.async_msgq = gr.msg_queue(0) - self.async_src = uhd.amsg_source("", self.async_msgq) - self.async_rcv = gru.msgq_runner(self.async_msgq, self.async_callback) - def async_callback(self, msg): - md = self.async_src.msg_to_async_metadata_t(msg) - print "Channel: %i Time: %f Event: %i" % (md.channel, md.time_spec.get_real_secs(), md.event_code) - self.error_code= md.event_code - self.stop() +# This function should behave exactly as MAIN, except it errors out +# as soon as any of the USRP errors are encountered. It should be run in +# a fashion like this: +# PYTHONPATH=. python -c "import peregrine.stream_usrp; peregrine.stream_usrp.main_capture_errors()" +# ... -1 -u name=MyB210 -d -g30 peregrine/sample_2015_09_11_18-47-11.1bit peregrine/a + +def main_capture_errors(): + args = sys.argv + args.pop(0) + args.insert(0,"peregrine/stream_usrp.py") + print args + proc = subprocess.Popen(args, + stderr=subprocess.PIPE) + out_str = "" + while proc.poll() == None: + errchar = proc.stderr.read(1) + if errchar == 'U': + print "Stream_usrp exiting due to Underflow at time {0}".format(str(datetime.datetime.now())) + proc.kill() + sys.exit(2) + if errchar == 'L': + print "Stream_usrp exiting due to Undeflow at time {0}".format(str(dateime.datetime.now())) + proc.kill() + sys.exit(3) + if errchar == "\n": + sys.stderr.write(out_str) + out_str = "" + else: + out_str += errchar + # Sleep for a second before exiting if it's not one of the cases we handle specially + time.sleep(1) + out_str += proc.stderr.read() + if out_str != "": + sys.stderr.write(out_str) + return proc.returncode def main(): if gr.enable_realtime_scheduling() != gr.RT_OK: @@ -179,7 +206,7 @@ def main(): help="Center frequency (%(default).0f)") parser.add_argument("-d", dest="dual", action='store_true', help="Using dual USRP devices") - parser.add_argument("-o", dest="outfile", default=None, + parser.add_argument("-o", dest="outfile", default="out.txt", help="Route Python stdout/stderr to this file") args = parser.parse_args() @@ -210,9 +237,7 @@ def main(): sync_pps=args.pps) tb.start() tb.wait() - if args.outfile: - stdout.close() - sys.exit(tb.error_code) + stdout.close() if __name__ == '__main__': main() diff --git a/peregrine/tracking.py b/peregrine/tracking.py index a94c216..15b88a8 100644 --- a/peregrine/tracking.py +++ b/peregrine/tracking.py @@ -1,4 +1,5 @@ -# Copyright (C) 2012 Swift Navigation Inc. +# Copyright (C) 2012,2016 Swift Navigation Inc. +# Contact: Adel Mamin # # This source is subject to the license found in the file 'LICENSE' which must # be be distributed together with this source. All other rights reserved. @@ -8,18 +9,30 @@ # WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A PARTICULAR PURPOSE. import numpy as np -from include.generateCAcode import caCodes -import gps_constants -import progressbar +import os import math import parallel_processing as pp - -import swiftnav.track -import swiftnav.correlate -import swiftnav.nav_msg -import defaults +import multiprocessing as mp +import cPickle + +from swiftnav.track import LockDetector +from swiftnav.track import CN0Estimator +from swiftnav.track import AliasDetector +from swiftnav.track import AidedTrackingLoop +from swiftnav.correlate import track_correlate +from swiftnav.nav_msg import NavMsg +from swiftnav.cnav_msg import CNavMsg +from swiftnav.cnav_msg import CNavMsgDecoder +from swiftnav.ephemeris import Ephemeris +from peregrine import defaults +from peregrine import gps_constants +from peregrine.acquisition import AcquisitionResult +from peregrine.include.generateCAcode import caCodes +from peregrine.include.generateL2CMcode import L2CMCodes import logging +import sys + logger = logging.getLogger(__name__) # Import progressbar if it is available. @@ -81,222 +94,631 @@ def update(self, e, p, l): """ raise NotImplementedError() -def track(samples, channels, - ms_to_track=None, - sampling_freq=defaults.sampling_freq, - chipping_rate=defaults.chipping_rate, - IF=defaults.IF, - show_progress=True, - loop_filter_class=swiftnav.track.AidedTrackingLoop, - stage1_loop_filter_params=( - (1, 0.7, 1), # Code loop NBW, zeta, k - (25, 0.7, 1), # Carrier loop NBW, zeta, k - 1e3, # Loop frequency - 5, # Carrier loop aiding_igain - 1540 - ), - correlator=swiftnav.correlate.track_correlate_, - stage2_coherent_ms=None, - stage2_loop_filter_params=None, - multi=True): - - n_channels = len(channels) - - # Add 22ms for safety, the corellator might try to access data a bit past - # just the number of milliseconds specified. - # TODO: Fix the correlator so this isn't an issue. - samples_length_ms = int(1e3 * len(samples) / sampling_freq - 22) - - if ms_to_track is None: - ms_to_track = samples_length_ms - - if samples_length_ms < ms_to_track: - logger.warning("Samples set too short for requested tracking length (%.4fs)" - % (ms_to_track * 1e-3)) - ms_to_track = samples_length_ms - - logger.info("Tracking %.4fs of data (%d samples)" % - (ms_to_track * 1e-3, ms_to_track * 1e-3 * sampling_freq)) - - # Make sure we have an integer number of points - num_points = int(math.floor(ms_to_track)) - - logger.info("Tracking starting") - logger.debug("Tracking %d channels, PRNs %s" % - (n_channels, [chan.prn+1 for chan in channels])) - - # If progressbar is not available, disable show_progress. - if show_progress and not _progressbar_available: - show_progress = False - logger.warning("show_progress = True but progressbar module not found.") - - # Setup our progress bar if we need it - if show_progress and not multi: - widgets = [' Tracking ', - progressbar.Attribute(['chan', 'nchan'], - '(CH: %d/%d)', - '(CH: -/-)'), ' ', - progressbar.Percentage(), ' ', - progressbar.ETA(), ' ', - progressbar.Bar()] - pbar = progressbar.ProgressBar(widgets=widgets, - maxval=n_channels*num_points, - attr={'nchan': n_channels}) - pbar.start() - else: - pbar = None - - # Run tracking for each channel - def do_channel(chan, n=None, q_progress=None): - track_result = TrackResults(num_points, chan.prn) - # Convert acquisition SNR to C/N0 - cn0_0 = 10 * np.log10(chan.snr) - cn0_0 += 10 * np.log10(1000) # Channel bandwidth - cn0_est = swiftnav.track.CN0Estimator( - bw=1e3, - cn0_0=cn0_0, - cutoff_freq=10, - loop_freq=1e3 - ) - - # Estimate initial code freq via aiding from acq carrier freq - code_freq_init = (chan.carr_freq - IF) * \ - gps_constants.chip_rate / gps_constants.l1 - carr_freq_init = chan.carr_freq - IF - loop_filter = loop_filter_class( - loop_freq = stage1_loop_filter_params[2], - code_freq = code_freq_init, - code_bw = stage1_loop_filter_params[0][0], - code_zeta = stage1_loop_filter_params[0][1], - code_k = stage1_loop_filter_params[0][2], - carr_to_code = stage1_loop_filter_params[4], - carr_freq = carr_freq_init, - carr_bw = stage1_loop_filter_params[1][0], - carr_zeta = stage1_loop_filter_params[1][1], - carr_k = stage1_loop_filter_params[1][2], - carr_freq_b1 = stage1_loop_filter_params[3], +def _tracking_channel_factory(parameters): + if parameters['acq'].signal == gps_constants.L1CA: + return TrackingChannelL1CA(parameters) + if parameters['acq'].signal == gps_constants.L2C: + return TrackingChannelL2C(parameters) + + +class TrackingChannel(object): + + def __init__(self, params): + for (key, value) in params.iteritems(): + setattr(self, key, value) + + self.prn = params['acq'].prn + self.signal = params['acq'].signal + + self.results_num = 500 + self.stage1 = True + + self.lock_detect = LockDetector( + k1=self.lock_detect_params["k1"], + k2=self.lock_detect_params["k2"], + lp=self.lock_detect_params["lp"], + lo=self.lock_detect_params["lo"]) + + self.alias_detect = AliasDetector( + acc_len=defaults.alias_detect_interval_ms / self.coherent_ms, + time_diff=1) + + self.cn0_est = CN0Estimator( + bw=1e3 / self.coherent_ms, + cn0_0=self.cn0_0, + cutoff_freq=10, + loop_freq=self.loop_filter_params["loop_freq"] ) - code_phase = 0.0 - carr_phase = 0.0 - # Get a vector with the C/A code sampled 1x/chip - ca_code = caCodes[chan.prn] + self.loop_filter = self.loop_filter_class( + loop_freq=self.loop_filter_params['loop_freq'], + code_freq=self.code_freq_init, + code_bw=self.loop_filter_params['code_bw'], + code_zeta=self.loop_filter_params['code_zeta'], + code_k=self.loop_filter_params['code_k'], + carr_to_code=0, # the provided code frequency accounts for Doppler + carr_freq=self.acq.doppler, + carr_bw=self.loop_filter_params['carr_bw'], + carr_zeta=self.loop_filter_params['carr_zeta'], + carr_k=self.loop_filter_params['carr_k'], + carr_freq_b1=self.loop_filter_params['carr_freq_b1'], + ) - # Add wrapping to either end to be able to do early/late - ca_code = np.concatenate(([ca_code[1022]], ca_code, [ca_code[0]])) + self.next_code_freq = self.loop_filter.to_dict()['code_freq'] + self.next_carr_freq = self.loop_filter.to_dict()['carr_freq'] + + self.track_result = TrackResults(self.results_num, + self.acq.prn, + self.acq.signal) + self.alias_detect_init = 1 + self.code_phase = 0.0 + self.carr_phase = 0.0 + self.samples_per_chip = int(round(self.sampling_freq / self.chipping_rate)) + self.sample_index = self.acq.sample_index + self.sample_index += self.acq.code_phase * self.samples_per_chip + self.sample_index = int(math.floor(self.sample_index)) + self.carr_phase_acc = 0.0 + self.code_phase_acc = 0.0 + self.samples_tracked = 0 + self.i = 0 + #self.samples_offset = self.samples['samples_offset'] + + self.pipelining = False # Flag if pipelining is used + self.pipelining_k = 0. # Error prediction coefficient for pipelining + self.short_n_long = False # Short/Long cycle simulation + self.short_step = False # Short cycle + if self.tracker_options: + mode = self.tracker_options['mode'] + if mode == 'pipelining': + self.pipelining = True + self.pipelining_k = self.tracker_options['k'] + elif mode == 'short-long-cycles': + self.short_n_long = True + self.pipelining = True + self.pipelining_k = self.tracker_options['k'] + else: + raise ValueError("Invalid tracker mode %s" % str(mode)) + + def dump(self): + fn_analysis, fn_results = self.track_result.dump(self.output_file, self.i) + self.i = 0 + return fn_analysis, fn_results + + def start(self): + logger.info("[PRN: %d (%s)] Tracking is started. " + "IF: %.1f, Doppler: %.1f, code phase: %.1f, " + "sample channel: %d sample index: %d" % + (self.prn + 1, + self.signal, + self.IF, + self.acq.doppler, + self.acq.code_phase, + self.acq.sample_channel, + self.acq.sample_index)) + + def get_index(self): + return self.sample_index + + def _run_preprocess(self): # optionally redefined in subclasses + pass + + def _run_postprocess(self): # optionally redefine in subclasses + pass + + def _get_result(self): # optionally redefine in subclasses + return None + + def is_pickleable(self): + return True - # Number of samples to seek ahead in file - samples_per_chip = int(round(sampling_freq / chipping_rate)) + def run_parallel(self, samples): + handover = self.run(samples) + return (handover, self) - # Set sample_index to start on a code rollover - sample_index = chan.code_phase * samples_per_chip + def run(self, samples): - # Start in 1ms integration until we know the nav bit phase - stage1 = True + self.samples = samples - carr_phase_acc = 0.0 - code_phase_acc = 0.0 + if self.sample_index < samples['sample_index']: + raise ValueError("Incorrent samples offset") - progress = 0 - ms_tracked = 0 - i = 0 - # Process the specified number of ms - while ms_tracked < ms_to_track: - if pbar: - pbar.update(ms_tracked + n * num_points, attr={'chan': n+1}) + sample_index = self.sample_index - samples['sample_index'] + samples_processed = 0 + samples_total = len(samples[self.signal]['samples']) - E = 0+0.j; P = 0+0.j; L = 0+0.j + estimated_blksize = self.coherent_ms * self.sampling_freq / 1e3 - if stage1 and stage2_coherent_ms and track_result.nav_msg.bit_phase == track_result.nav_msg.bit_phase_ref: - #print "PRN %02d transition to stage 2 at %d ms" % (chan.prn+1, ms_tracked) - stage1 = False - loop_filter.retune(*stage2_loop_filter_params) - cn0_est = swiftnav.track.CN0Estimator(bw=1e3/stage2_coherent_ms, - cn0_0=track_result.cn0[i-1], - cutoff_freq=10, - loop_freq=1e3/stage2_coherent_ms) + while self.samples_tracked < self.samples_to_track and \ + (sample_index + 2 * estimated_blksize) < samples_total: - coherent_ms = 1 if stage1 else stage2_coherent_ms + self._run_preprocess() - for j in range(coherent_ms): - samples_ = samples[sample_index:] + if self.pipelining: + # Pipelining and prediction + corr_code_freq, corr_carr_freq = self.next_code_freq, self.next_carr_freq - E_, P_, L_, blksize, code_phase, carr_phase = correlator( - samples_, - loop_filter.to_dict()['code_freq'] + chipping_rate, code_phase, - loop_filter.to_dict()['carr_freq'] + IF, carr_phase, - ca_code, - sampling_freq - ) - sample_index += blksize - carr_phase_acc += loop_filter.to_dict()['carr_freq'] * blksize / sampling_freq - code_phase_acc += loop_filter.to_dict()['code_freq'] * blksize / sampling_freq + self.next_code_freq = self.loop_filter.to_dict()['code_freq'] + self.next_carr_freq = self.loop_filter.to_dict()['carr_freq'] + + if self.short_n_long and not self.stage1 and not self.short_step: + # In case of short/long cycles, the correction applicable for the + # long cycle is smaller proportionally to the actual cycle size + pipelining_k = self.pipelining_k / (self.coherent_ms - 1) + else: + pipelining_k = self.pipelining_k + + # There is an error between target frequency and actual one. Affect + # the target frequency according to the computed error + carr_freq_error = self.next_carr_freq - corr_carr_freq + self.next_carr_freq += carr_freq_error * pipelining_k + + code_freq_error = self.next_code_freq - corr_code_freq + self.next_code_freq += code_freq_error * pipelining_k - E += E_; P += P_; L += L_ + else: + # Immediate correction simulation + self.next_code_freq = self.loop_filter.to_dict()['code_freq'] + self.next_carr_freq = self.loop_filter.to_dict()['carr_freq'] - loop_filter.update(E, P, L) - track_result.coherent_ms[i] = coherent_ms + corr_code_freq, corr_carr_freq = self.next_code_freq, self.next_carr_freq - track_result.nav_bit_sync.update(np.real(P), coherent_ms) + if self.short_n_long and not self.stage1: + # When simulating short and long cycles, short step resets EPL + # registers, and long one adds up to them + if self.short_step: + self.E = self.P = self.L = 0.j + self.coherent_iter = 1 + else: + self.coherent_iter = self.coherent_ms - 1 + else: + self.E = self.P = self.L = 0.j - # TODO - Is this the correct way to call nav_msg.update? - tow = track_result.nav_msg.update(np.real(P) >= 0) - track_result.nav_msg_bit_phase_ref[i] = track_result.nav_msg.bit_phase_ref - track_result.tow[i] = tow or (track_result.tow[i-1] + coherent_ms) + for _ in range(self.coherent_iter): - track_result.carr_phase[i] = carr_phase - track_result.carr_phase_acc[i] = carr_phase_acc - track_result.carr_freq[i] = loop_filter.to_dict()['carr_freq'] + IF + if (sample_index + 2 * estimated_blksize) >= samples_total: + break - track_result.code_phase[i] = code_phase - track_result.code_phase_acc[i] = code_phase_acc - track_result.code_freq[i] = loop_filter.to_dict()['code_freq'] + chipping_rate + samples_ = samples[self.signal]['samples'][sample_index:] + + E_, P_, L_, blksize, self.code_phase, self.carr_phase = self.correlator( + samples_, + corr_code_freq + self.chipping_rate, self.code_phase, + corr_carr_freq + self.IF, self.carr_phase, + self.prn_code, + self.sampling_freq, + self.signal + ) + + if blksize > estimated_blksize: + estimated_blksize = blksize + + sample_index += blksize + samples_processed += blksize + self.carr_phase_acc += corr_carr_freq * blksize / self.sampling_freq + self.code_phase_acc += corr_code_freq * blksize / self.sampling_freq + + self.E += E_ + self.P += P_ + self.L += L_ + + if not self.stage1 and self.short_n_long: + if self.short_step: + # In case of short step - go to next integration period + self.short_step = False + self.alias_detect.first(self.P.real, self.P.imag) + continue + else: + # Next step is short cycle + self.short_step = True + + # Update PLL lock detector + lock_detect_outo, \ + lock_detect_outp, \ + lock_detect_pcount1, \ + lock_detect_pcount2, \ + lock_detect_lpfi, \ + lock_detect_lpfq = self.lock_detect.update(self.P.real, + self.P.imag, + self.coherent_iter) + + if lock_detect_outo: + if self.alias_detect_init: + self.alias_detect_init = 0 + self.alias_detect.reinit(defaults.alias_detect_interval_ms / + self.coherent_iter, + time_diff=1) + self.alias_detect.first(self.P.real, self.P.imag) + alias_detect_err_hz = \ + self.alias_detect.second(self.P.real, self.P.imag) * np.pi * \ + (1e3 / defaults.alias_detect_interval_ms) + self.alias_detect.first(self.P.real, self.P.imag) + else: + self.alias_detect_init = 1 + alias_detect_err_hz = 0 + + self.loop_filter.update(self.E, self.P, self.L) + self.track_result.coherent_ms[self.i] = self.coherent_ms + + self.track_result.IF = self.IF + self.track_result.carr_phase[self.i] = self.carr_phase + self.track_result.carr_phase_acc[self.i] = self.carr_phase_acc + self.track_result.carr_freq[self.i] = \ + self.loop_filter.to_dict()['carr_freq'] + self.IF + + self.track_result.code_phase[self.i] = self.code_phase + self.track_result.code_phase_acc[self.i] = self.code_phase_acc + self.track_result.code_freq[self.i] = \ + self.loop_filter.to_dict()['code_freq'] + self.chipping_rate # Record stuff for postprocessing - track_result.absolute_sample[i] = sample_index + self.track_result.absolute_sample[self.i] = self.sample_index + \ + samples_processed - track_result.E[i] = E - track_result.P[i] = P - track_result.L[i] = L + self.track_result.E[self.i] = self.E + self.track_result.P[self.i] = self.P + self.track_result.L[self.i] = self.L - track_result.cn0[i] = cn0_est.update(P.real, P.imag) + self.track_result.cn0[self.i] = self.cn0_est.update( + self.P.real, self.P.imag) - i += 1 - ms_tracked += coherent_ms + self.track_result.lock_detect_outo[self.i] = lock_detect_outo + self.track_result.lock_detect_outp[self.i] = lock_detect_outp + self.track_result.lock_detect_pcount1[self.i] = lock_detect_pcount1 + self.track_result.lock_detect_pcount2[self.i] = lock_detect_pcount2 + self.track_result.lock_detect_lpfi[self.i] = lock_detect_lpfi + self.track_result.lock_detect_lpfq[self.i] = lock_detect_lpfq - if q_progress and (i % 200 == 0): - p = 1.0 * ms_tracked / ms_to_track - q_progress.put(p - progress) - progress = p + self.track_result.alias_detect_err_hz[self.i] = alias_detect_err_hz - # Possibility for lock-detection later - track_result.status = 'T' + self._run_postprocess() - track_result.resize(i) - if q_progress: - q_progress.put(1.0 - progress) + self.samples_tracked = self.sample_index + samples_processed + self.track_result.ms_tracked[self.i] = self.samples_tracked * 1e3 / \ + self.sampling_freq - return track_result + self.i += 1 + if self.i >= self.results_num: + self.dump() - if multi: - track_results=pp.parmap(do_channel, channels, - show_progress=show_progress, func_progress=show_progress) - else: - track_results=map(lambda (n, chan): do_channel(chan, n=n), enumerate(channels)) + self.sample_index += samples_processed + self.track_result.status = 'T' - if pbar: - pbar.finish() + return self._get_result() - logger.info("Tracking finished") - return track_results +class TrackingChannelL1CA(TrackingChannel): + + def __init__(self, params): + # Convert acquisition SNR to C/N0 + cn0_0 = 10 * np.log10(params['acq'].snr) + cn0_0 += 10 * np.log10(defaults.L1CA_CHANNEL_BANDWIDTH_HZ) + + params['cn0_0'] = cn0_0 + params['coherent_ms'] = 1 + params['IF'] = params['samples'][gps_constants.L1CA]['IF'] + params['prn_code'] = caCodes[params['acq'].prn] + params['code_freq_init'] = params['acq'].doppler * \ + gps_constants.chip_rate / gps_constants.l1 + params['loop_filter_params'] = defaults.l1ca_stage1_loop_filter_params + params['lock_detect_params'] = defaults.l1ca_lock_detect_params_opt + + TrackingChannel.__init__(self, params) + + self.nav_msg = NavMsg() + self.nav_bit_sync = NBSMatchBit() if self.prn < 32 else NBSSBAS() + self.l2c_handover_acq = None + self.l2c_handover_done = False + + def _run_preprocess(self): + # For L1 C/A there are coherent and non-coherent tracking options. + if self.stage1 and \ + self.stage2_coherent_ms and \ + self.nav_bit_sync.bit_phase == self.nav_bit_sync.bit_phase_ref: + + logger.info("[PRN: %d (%s)] switching to stage2, coherent_ms=%d" % + (self.prn + 1, self.signal, self.stage2_coherent_ms)) + + self.stage1 = False + self.coherent_ms = self.stage2_coherent_ms + + self.loop_filter.retune(**self.stage2_loop_filter_params) + self.lock_detect.reinit( + k1=self.lock_detect_params["k1"] * self.coherent_ms, + k2=self.lock_detect_params["k2"], + lp=self.lock_detect_params["lp"], + lo=self.lock_detect_params["lo"]) + self.cn0_est = CN0Estimator(bw=1e3 / self.stage2_coherent_ms, + cn0_0=self.track_result.cn0[self.i - 1], + cutoff_freq=10, + loop_freq=1e3 / self.stage2_coherent_ms) + + self.coherent_iter = self.coherent_ms + + def _get_result(self): + if self.l2c_handover_acq and not self.l2c_handover_done: + self.l2c_handover_done = True + return self.l2c_handover_acq + return None + + def _run_postprocess(self): + sync, bit = self.nav_bit_sync.update(np.real(self.P), self.coherent_ms) + if sync: + tow = self.nav_msg.update(bit) + if tow >= 0: + logger.info("[PRN: %d (%s)] ToW %d" % + (self.prn + 1, self.signal, tow)) + if self.nav_msg.subframe_ready(): + eph = Ephemeris() + res = self.nav_msg.process_subframe(eph) + if res < 0: + logger.error("[PRN: %d (%s)] Subframe decoding error %d" % + (self.prn + 1, self.signal, res)) + elif res > 0: + logger.info("[PRN: %d (%s)] Subframe decoded" % + (self.prn + 1, self.signal)) + else: + # Subframe decoding is in progress + pass + else: + tow = -1 + self.track_result.tow[self.i] = tow if tow >= 0 else ( + self.track_result.tow[self.i - 1] + self.coherent_ms) + + # Handover to L2C if possible + if self.l2c_handover and not self.l2c_handover_acq and \ + 'samples' in self.samples[gps_constants.L2C] and sync: + chan_snr = self.track_result.cn0[self.i] + chan_snr -= 10 * np.log10(defaults.L1CA_CHANNEL_BANDWIDTH_HZ) + chan_snr = np.power(10, chan_snr / 10) + l2c_doppler = self.loop_filter.to_dict( + )['carr_freq'] * gps_constants.l2 / gps_constants.l1 + self.l2c_handover_acq = \ + AcquisitionResult(self.prn, + self.samples[gps_constants.L2C][ + 'IF'] + l2c_doppler, + l2c_doppler, # carrier doppler + self.track_result.code_phase[ + self.i], + chan_snr, + 'A', + gps_constants.L2C, + self.track_result.absolute_sample[self.i]) + + +class TrackingChannelL2C(TrackingChannel): + + def __init__(self, params): + # Convert acquisition SNR to C/N0 + cn0_0 = 10 * np.log10(params['acq'].snr) + cn0_0 += 10 * np.log10(defaults.L2C_CHANNEL_BANDWIDTH_HZ) + params['cn0_0'] = cn0_0 + params['coherent_ms'] = 20 + params['coherent_iter'] = 1 + params['loop_filter_params'] = defaults.l2c_loop_filter_params + params['lock_detect_params'] = defaults.l2c_lock_detect_params_20ms + params['IF'] = params['samples'][gps_constants.L2C]['IF'] + params['prn_code'] = L2CMCodes[params['acq'].prn] + params['code_freq_init'] = params['acq'].doppler * \ + gps_constants.chip_rate / gps_constants.l2 + + TrackingChannel.__init__(self, params) + + self.cnav_msg = CNavMsg() + self.cnav_msg_decoder = CNavMsgDecoder() + + def is_pickleable(self): + return False # could not pickle cnav_msg_decoder Cython object + + def _run_postprocess(self): + symbol = 0xFF if np.real(self.P) >= 0 else 0x00 + res, delay = self.cnav_msg_decoder.decode(symbol, self.cnav_msg) + if res: + logger.debug("[PRN: %d (%s)] CNAV message decoded: " + "prn=%d msg_id=%d tow=%d alert=%d delay=%d" % + (self.prn + 1, + self.signal, + self.cnav_msg.getPrn(), + self.cnav_msg.getMsgId(), + self.cnav_msg.getTow(), + self.cnav_msg.getAlert(), + delay)) + tow = self.cnav_msg.getTow() * 6000 + delay * 20 + logger.debug("[PRN: %d (%s)] ToW %d" % + (self.prn + 1, self.signal, tow)) + self.track_result.tow[self.i] = tow + else: + self.track_result.tow[self.i] = self.track_result.tow[self.i - 1] + \ + self.coherent_ms + + +class Tracker(object): + + def __init__(self, + samples, + channels, + ms_to_track, + sampling_freq, + chipping_rate=defaults.chipping_rate, + l2c_handover=True, + progress_bar_output='none', + loop_filter_class=AidedTrackingLoop, + correlator=track_correlate, + stage2_coherent_ms=None, + stage2_loop_filter_params=None, + multi=False, + tracker_options=None, + output_file=None): + + self.samples = samples + self.sampling_freq = sampling_freq + self.ms_to_track = ms_to_track + self.tracker_options = tracker_options + self.output_file = output_file + self.chipping_rate = chipping_rate + self.l2c_handover = l2c_handover + self.correlator = correlator + self.stage2_coherent_ms = stage2_coherent_ms + self.stage2_loop_filter_params = stage2_loop_filter_params + + if mp.cpu_count() > 1: + self.multi = multi + else: + self.multi = False + + self.loop_filter_class = loop_filter_class + + if self.ms_to_track: + self.samples_to_track = self.ms_to_track * sampling_freq / 1e3 + if samples['samples_total'] < self.samples_to_track: + logger.warning("Samples set too short for requested tracking length (%.4fs)" + % (self.ms_to_track * 1e-3)) + self.samples_to_track = samples['samples_total'] + else: + self.samples_to_track = samples['samples_total'] + + if progress_bar_output == 'stdout': + self.show_progress = True + progress_fd = sys.stdout + elif progress_bar_output == 'stderr': + self.show_progress = True + progress_fd = sys.stderr + else: + self.show_progress = False + progress_fd = -1 + + # If progressbar is not available, disable show_progress. + if self.show_progress and not _progressbar_available: + self.show_progress = False + logger.warning("show_progress = True but progressbar module not found.") + + # Setup our progress bar if we need it + if self.show_progress: + widgets = [' Tracking ', + progressbar.Attribute(['sample', 'samples'], + '(sample: %d/%d)', + '(sample: -/-)'), ' ', + progressbar.Percentage(), ' ', + progressbar.ETA(), ' ', + progressbar.Bar()] + self.pbar = progressbar.ProgressBar(widgets=widgets, + maxval=samples['samples_total'], + attr={'samples': self.samples['samples_total'], + 'sample': 0l}, + fd=progress_fd) + else: + self.pbar = None + + self.tracking_channels = map(self._create_channel, channels) + + def start(self): + logger.info("Number of CPUs: %d" % (mp.cpu_count())) + + logger.info("Tracking %.4fs of data (%d samples)" % + (self.samples_to_track / self.sampling_freq, + self.samples_to_track)) + + logger.info("Tracking starting") + logger.debug("Tracking PRNs %s" % + ([chan.prn + 1 for chan in self.tracking_channels])) + + if self.pbar: + self.pbar.start() + + def _print_name(self, name): + print name + + def stop(self): + if self.pbar: + self.pbar.finish() + + # (fn_analysis, fn_results) = map(lambda chan: chan.dump(), self.tracking_channels) + res = map(lambda chan: chan.dump(), self.tracking_channels) + + fn_analysis = map(lambda x: x[0], res) + fn_results = map(lambda x: x[1], res) + + print "The tracking results were stored into:" + map(self._print_name, fn_analysis) + + logger.info("Tracking finished") + + return fn_results + + def _create_channel(self, acq): + if not acq: + return + parameters = {'acq': acq, + 'samples': self.samples, + 'loop_filter_class': self.loop_filter_class, + 'tracker_options': self.tracker_options, + 'output_file': self.output_file, + 'samples_to_track': self.samples_to_track, + 'sampling_freq': self.sampling_freq, + 'chipping_rate': self.chipping_rate, + 'l2c_handover': self.l2c_handover, + 'show_progress': self.show_progress, + 'correlator': self.correlator, + 'stage2_coherent_ms': self.stage2_coherent_ms, + 'stage2_loop_filter_params': self.stage2_loop_filter_params, + 'multi': self.multi} + return _tracking_channel_factory(parameters) + + def _run_parallel(self, i, samples): + handover = self.parallel_channels[i].run(samples) + return self.parallel_channels[i], handover + + def run_channels(self, samples): + channels = self.tracking_channels + self.tracking_channels = [] + + while channels and not all(v is None for v in channels): + + if self.multi: + self.parallel_channels = filter(lambda x: x.is_pickleable(), channels) + else: + self.parallel_channels = [] + + serial_channels = list(set(channels) - set(self.parallel_channels)) + channels = [] + handover = [] + + if self.parallel_channels: + res = pp.parmap(lambda i: self._run_parallel(i, samples), + range(len(self.parallel_channels)), + nprocs = len(self.parallel_channels), + show_progress=False, + func_progress=False) + + channels = map(lambda x: x[0], res) + handover += map(lambda x: x[1], res) + + if serial_channels: + handover += map(lambda x: x.run(samples), serial_channels) + + self.tracking_channels += channels + serial_channels + handover = [h for h in handover if h is not None] + if handover: + channels = map(self._create_channel, handover) + else: + channels = None + + indicies = map(lambda x: x.get_index(), self.tracking_channels) + min_index = min(indicies) + + if self.pbar: + self.pbar.update(min_index, attr={'sample': min_index}) + + return min_index class TrackResults: - def __init__(self, n_points, prn): + + def __init__(self, n_points, prn, signal): + self.print_start = 1 self.status = '-' self.prn = prn + self.IF = 0 self.absolute_sample = np.zeros(n_points) self.code_phase = np.zeros(n_points) self.code_phase_acc = np.zeros(n_points) @@ -308,12 +730,83 @@ def __init__(self, n_points, prn): self.P = np.zeros(n_points, dtype=np.complex128) self.L = np.zeros(n_points, dtype=np.complex128) self.cn0 = np.zeros(n_points) - self.nav_msg = swiftnav.nav_msg.NavMsg() + self.lock_detect_outp = np.zeros(n_points) + self.lock_detect_outo = np.zeros(n_points) + self.lock_detect_pcount1 = np.zeros(n_points) + self.lock_detect_pcount2 = np.zeros(n_points) + self.lock_detect_lpfi = np.zeros(n_points) + self.lock_detect_lpfq = np.zeros(n_points) + self.alias_detect_err_hz = np.zeros(n_points) + self.nav_msg = NavMsg() self.nav_msg_bit_phase_ref = np.zeros(n_points) self.nav_bit_sync = NBSMatchBit() if prn < 32 else NBSSBAS() self.tow = np.empty(n_points) self.tow[:] = np.NAN self.coherent_ms = np.zeros(n_points) + # self.cnav_msg = swiftnav.cnav_msg.CNavMsg() + # self.cnav_msg_decoder = swiftnav.cnav_msg.CNavMsgDecoder() + self.signal = signal + self.ms_tracked = np.zeros(n_points) + + def dump(self, output_file, size): + output_filename, output_file_extension = os.path.splitext(output_file) + + # mangle the analyses file name with the tracked signal name + fn_analysis = output_filename + \ + (".PRN-%d.%s" % (self.prn + 1, self.signal)) +\ + output_file_extension + + # mangle the results file name with the tracked signal name + fn_results = output_filename + \ + (".PRN-%d.%s" % (self.prn + 1, self.signal)) +\ + output_file_extension + '.track_results' + + if self.print_start: + mode = 'w' + else: + mode = 'a' + + # saving tracking results for navigation stage + with open(fn_results, mode) as f1: + cPickle.dump(self, f1, protocol = cPickle.HIGHEST_PROTOCOL) + + with open(fn_analysis, mode) as f1: + if self.print_start: + f1.write("sample_index,ms_tracked,IF,doppler_phase,carr_doppler," + "code_phase, code_freq," + "CN0,E_I,E_Q,P_I,P_Q,L_I,L_Q," + "lock_detect_outp,lock_detect_outo," + "lock_detect_pcount1,lock_detect_pcount2," + "lock_detect_lpfi,lock_detect_lpfq,alias_detect_err_hz," + "code_phase_acc\n") + for i in range(size): + f1.write("%s," % int(self.absolute_sample[i])) + f1.write("%s," % self.ms_tracked[i]) + f1.write("%s," % self.IF) + f1.write("%s," % self.carr_phase[i]) + f1.write("%s," % (self.carr_freq[i] - + self.IF)) + f1.write("%s," % self.code_phase[i]) + f1.write("%s," % self.code_freq[i]) + f1.write("%s," % self.cn0[i]) + f1.write("%s," % self.E[i].real) + f1.write("%s," % self.E[i].imag) + f1.write("%s," % self.P[i].real) + f1.write("%s," % self.P[i].imag) + f1.write("%s," % self.L[i].real) + f1.write("%s," % self.L[i].imag) + f1.write("%s," % int(self.lock_detect_outp[i])) + f1.write("%s," % int(self.lock_detect_outo[i])) + f1.write("%s," % int(self.lock_detect_pcount1[i])) + f1.write("%s," % int(self.lock_detect_pcount2[i])) + f1.write("%s," % self.lock_detect_lpfi[i]) + f1.write("%s," % self.lock_detect_lpfq[i]) + f1.write("%s," % self.alias_detect_err_hz[i]) + f1.write("%s\n" % self.code_phase_acc[i]) + + self.print_start = 0 + + return fn_analysis, fn_results def resize(self, n_points): for k in dir(self): @@ -341,11 +834,11 @@ def _equal(self, other): """ if self.__dict__.keys() != other.__dict__.keys(): return False - + for k in self.__dict__.keys(): if isinstance(self.__dict__[k], np.ndarray): # If np.ndarray, elements might be floats, so compare accordingly. - if any(np.greater((self.__dict__[k]-other.__dict__[k]), np.ones(len(self.__dict__[k]))*10e-6)): + if any(np.greater((self.__dict__[k] - other.__dict__[k]), np.ones(len(self.__dict__[k])) * 10e-6)): return False elif self.__dict__[k] != other.__dict__[k]: return False @@ -354,14 +847,15 @@ def _equal(self, other): class NavBitSync: + def __init__(self): self.bit_phase = 0 self.bit_integrate = 0 - self.synced=False - self.bits=[] - self.bit_phase_ref=-1 # A new bit begins when bit_phase == bit_phase_ref + self.synced = False + self.bits = [] + self.bit_phase_ref = -1 # A new bit begins when bit_phase == bit_phase_ref self.count = 0 - + def update(self, corr, ms): self.bit_phase += ms self.bit_phase %= 20 @@ -370,8 +864,12 @@ def update(self, corr, ms): if not self.synced: self.update_bit_sync(corr, ms) if self.bit_phase == self.bit_phase_ref: - self.bits.append(1 if self.bit_integrate > 0 else 0) + bit = 1 if self.bit_integrate > 0 else 0 + self.bits.append(bit) self.bit_integrate = 0 + return True, bit + else: + return False, None def update_bit_sync(self, corr, ms): raise NotImplementedError @@ -402,24 +900,26 @@ def _equal(self, other): """ if self.__dict__.keys() != other.__dict__.keys(): return False - + for k in self.__dict__.keys(): if isinstance(self.__dict__[k], np.ndarray): # If np.ndarray, elements might be floats, so compare accordingly. - if any((self.__dict__[k]-other.__dict__[k]) > 10e-6): + if any((self.__dict__[k] - other.__dict__[k]) > 10e-6): return False elif self.__dict__[k] != other.__dict__[k]: return False return True + class NavBitSyncSBAS: + def __init__(self): self.bit_phase = 0 self.bit_integrate = 0 - self.synced=False - self.bits=[] - self.bit_phase_ref=-1 # A new bit begins when bit_phase == bit_phase_ref + self.synced = False + self.bits = [] + self.bit_phase_ref = -1 # A new bit begins when bit_phase == bit_phase_ref self.count = 0 def update(self, corr, ms): @@ -441,6 +941,7 @@ def bitstring(self): class NBSSBAS(NavBitSyncSBAS): + def __init__(self, thres=200): NavBitSyncSBAS.__init__(self) self.hist = np.zeros(2) @@ -464,18 +965,22 @@ def update_bit_sync(self, corr, ms): self.synced = True self.bit_phase_ref = np.argmax(self.hist) + class NBSLibSwiftNav(NavBitSync): + def __init__(self): NavBitSync.__init__(self) - self.nav_msg = swiftnav.nav_msg.NavMsg() + self.nav_msg = NavMsg() def update_bit_sync(self, corr, ms): self.nav_msg.update(corr, ms) self.bit_phase_ref = self.nav_msg.bit_phase_ref self.synced = self.bit_phase_ref >= 0 + class NBSMatchBit(NavBitSync): - def __init__(self, thres=22): + + def __init__(self, thres=25): NavBitSync.__init__(self) self.hist = np.zeros(20) self.acc = 0 @@ -498,7 +1003,9 @@ def update_bit_sync(self, corr, ms): self.synced = True self.bit_phase_ref = np.argmax(self.hist) + class NBSHistogram(NavBitSync): + def __init__(self, thres=10): NavBitSync.__init__(self) self.bit_phase_count = 0 @@ -518,8 +1025,11 @@ def update_bit_sync(self, corr, ms): self.hist = np.zeros(20) self.bit_phase_count = 0 + class NBSMatchEdge(NavBitSync): - # TODO: This isn't quite right - might get wrong answer with long leading run of same bits, depending on initial phase + # TODO: This isn't quite right - might get wrong answer with long leading + # run of same bits, depending on initial phase + def __init__(self, thres=100000): NavBitSync.__init__(self) self.hist = np.zeros(20) @@ -529,7 +1039,7 @@ def __init__(self, thres=100000): def update_bit_sync(self, corr, ms): bp40 = self.bit_phase % 40 - self.acc += corr - 2*self.prev[(bp40 - 20) % 40] + self.prev[bp40] + self.acc += corr - 2 * self.prev[(bp40 - 20) % 40] + self.prev[bp40] self.prev[bp40] = corr if self.bit_phase >= 40: # Accumulator valid diff --git a/tests/gpsl1ca_ci_samples.piksi_format.acq_results b/tests/gpsl1ca_ci_samples.piksi_format.acq_results new file mode 100644 index 0000000..920e154 Binary files /dev/null and b/tests/gpsl1ca_ci_samples.piksi_format.acq_results differ diff --git a/tests/test_generateL2CMcode.py b/tests/test_generateL2CMcode.py new file mode 100644 index 0000000..b682449 --- /dev/null +++ b/tests/test_generateL2CMcode.py @@ -0,0 +1,42 @@ +import pytest +from peregrine.include.generateL2CMcode import generateL2CMcode + +end_shift_regs_test = [\ + 0552566002, + 0034445034, + 0723443711, + 0511222013, + 0463055213, + 0667044524, + 0652322653, + 0505703344, + 0520302775, + 0244205506, + 0236174002, + 0654305531, + 0435070571, + 0630431251, + 0234043417, + 0535540745, + 0043056734, + 0731304103, + 0412120105, + 0365636111, + 0143324657, + 0110766462, + 0602405203, + 0177735650, + 0630177560, + 0653467107, + 0406576630, + 0221777100, + 0773266673, + 0100010710, + 0431037132, + 0624127475 +] + +def test_generateL2CMcode(): + for i in range(32): + assert (end_shift_regs_test[i] == generateL2CMcode(i)[1]) + diff --git a/tests/test_run.py b/tests/test_run.py index f7cf3b3..17c0783 100644 --- a/tests/test_run.py +++ b/tests/test_run.py @@ -17,7 +17,12 @@ from shutil import copyfile SAMPLES_PATH = 'tests/test_data/' -RES_PATH = SAMPLES_PATH + '/results/' +# todo: the gpsl1ca_ci_samples.piksi_format.acq_results +# should replace the old file with the same name at the +# remote server, where the script takes it from. +# For now, let's use the local version. +#RES_PATH = SAMPLES_PATH + '/results/' +RES_PATH = 'tests/' SAMPLES_FNAME = 'gpsl1ca_ci_samples.piksi_format' @@ -77,5 +82,5 @@ def test_tracking(): # Clean-up. os.remove(NEW_ACQ_RES) - os.remove(NEW_TRK_RES) + #os.remove(NEW_TRK_RES)