From 90d50c53c35fa615dc7e6ff180d038d46e06678b Mon Sep 17 00:00:00 2001 From: Perttu Salmela Date: Sun, 5 Jun 2016 20:44:33 +0300 Subject: [PATCH] Extend performance simulation scripts with L1CA band option --- peregrine/analysis/README-sim.txt | 20 ++- peregrine/analysis/plt_res.py | 56 +++++-- peregrine/analysis/run_sim.py | 233 +++++++++++++++++++++--------- 3 files changed, 220 insertions(+), 89 deletions(-) diff --git a/peregrine/analysis/README-sim.txt b/peregrine/analysis/README-sim.txt index c210db4..d62f87c 100644 --- a/peregrine/analysis/README-sim.txt +++ b/peregrine/analysis/README-sim.txt @@ -1,14 +1,14 @@ -Running L2C Performance Simulations +Running L1CA/L2C Performance Simulations =================================== Performance simulations are run with run_sim.py tool which invokes iqgen_main.py and tracking_loop.py. The tool stores all the simulation results (including commands -for generating iq data and tracking l2c signal) in json file. Data is always appended +for generating iq data and tracking l1ca/l2c signal) in json file. Data is always appended in given json file. Thereafter, results are illustrated with plt_res.py which reads data from json file. -The json file contains following data for each iqgen and L2C tracking run: +The json file contains following data for each iqgen and L1CA/L2C tracking run: { "acc": # Acceleration (L1CA Hz / s) "avgcn0": # Average of all tracker CN0 estimates @@ -18,9 +18,9 @@ The json file contains following data for each iqgen and L2C tracking run: "dopSigma3": # 3-sigma Doppler error (99.7% of errors is less than this) "duration": # Length of simulation in seconds "iqgencmd": # Command line for iqgen - "iqgencn0": # L2C CN0 reported by iqgen - "l2chip": # Initial L2C code phase reported by iqgen - "l2dop": # Initial L2C Doppler frequency reported by iqgen + "iqgencn0": # L1CA/L2C CN0 reported by iqgen + "l_chip": # Initial L1CA/L2C code phase reported by iqgen + "l_dop": # Initial L1CA/L2C Doppler frequency reported by iqgen "lockrate": # PLL lock rate "snr": # SNR argument for iqgen "stamp": # Wall clock time stamp of simulation run @@ -52,6 +52,13 @@ Performance data can be generated with following commands: $ ./plt_res.py -e -f result_3.json +Band selection +============== +By default simulation is run for L2C band. The L1CA is selected with command line option: +-b l1ca (--band) +FPGA delay control simulation option --short-long-cycles can be given for L1CA band simulation with option: +-s (--short-long-cycles) + Varying simulation parameters ============================= Changing values in the python file is as easy as value in configuration file @@ -89,3 +96,4 @@ Changing CN0 estimator of plt_res.py In plt_res.py there is a string which defines which CN0 estimate is used: CN0STRING="avgcn0" # CN0 from tracking #CN0STRING="iqgencn0" # CN0 from iqgen + diff --git a/peregrine/analysis/plt_res.py b/peregrine/analysis/plt_res.py index cc47c06..c4c4c47 100755 --- a/peregrine/analysis/plt_res.py +++ b/peregrine/analysis/plt_res.py @@ -16,9 +16,22 @@ import json import numpy import matplotlib.pyplot as plt +import peregrine.gps_constants -CN0STRING = "avgcn0" # CN0 from tracking -# CN0STRING="iqgencn0" # CN0 from iqgen + +class CfgClass: + + def __init__(self): + # self.CN0STRING = "avgcn0" # CN0 from tracking + self.CN0STRING = "iqgencn0" # CN0 from iqgen + self.BAND = "l2c" + + def isL1CA(self): + if "l1ca" == self.BAND: + return True + return False + +cfg = CfgClass() def sigmaFreqPlot(filename): @@ -30,7 +43,7 @@ def sigmaFreqPlot(filename): jsarray = json.loads(s) r = [] for j in jsarray: - r.append((j["dopSigma1"], j[CN0STRING])) + r.append((j["dopSigma1"], j[cfg.CN0STRING])) r = sorted(r, key=lambda x: x[1]) dopSigma1 = map(lambda x: x[0], r) @@ -40,6 +53,7 @@ def sigmaFreqPlot(filename): plt.plot(avgcn0, dopSigma1, 'o-') plt.xlabel('CN0') plt.ylabel('Doppler sigma-1 error (Hz)') + plt.grid(True) plt.show() @@ -52,7 +66,7 @@ def lockrateCn0Plot(filename): jsarray = json.loads(s) r = [] for j in jsarray: - r.append((j["lockrate"], j[CN0STRING])) + r.append((j["lockrate"], j[cfg.CN0STRING])) r = sorted(r, key=lambda x: x[1]) lockrate = map(lambda x: x[0], r) @@ -63,10 +77,14 @@ def lockrateCn0Plot(filename): plt.xlabel('CN0') plt.ylabel('PLL lock rate') #plt.scatter(avgcn0, lockrate) + plt.grid(True) plt.show() def dynamicAccPlot(filename, mode): + GRAV_G = 9.80665 # m/s^2 + HzToMps = (peregrine.gps_constants.c / peregrine.gps_constants.l1) + HzToG = HzToMps / GRAV_G fp = open(filename, "r") s = fp.read() fp.close() @@ -78,8 +96,8 @@ def dynamicAccPlot(filename, mode): maxCN0 = 0.0 r = [] for j in jsarray: - minCN0 = min(minCN0, j[CN0STRING]) - maxCN0 = max(maxCN0, j[CN0STRING]) + minCN0 = min(minCN0, j[cfg.CN0STRING]) + maxCN0 = max(maxCN0, j[cfg.CN0STRING]) #r.append( (float(j["acc"]),j[CN0STRING]) ) #accVecAll = map(lambda x:x[0], r) #cn0VecAll = map(lambda x:x[1], r) @@ -89,23 +107,28 @@ def dynamicAccPlot(filename, mode): fig = plt.figure() r = [] + if cfg.isL1CA(): + Tcoh = 0.005 # Integration time 5 ms + else: + Tcoh = 0.02 # Integration time 20 ms + for cn0bin in cn0Range: bestAcc = -1000.0 bestCN0 = 0 print "BIN", cn0bin, for j in jsarray: - cn0 = j[CN0STRING] + cn0 = j[cfg.CN0STRING] lockrate = j["lockrate"] doperr = j["dopSigma1"] - acc = float(j["acc"]) + acc = float(j["acc"]) * HzToG if cn0bin - 0.5 <= cn0 and cn0 < cn0bin + 0.5: if acc > bestAcc: if mode == "lockrate" and lockrate >= 0.68: bestAcc = acc bestCN0 = cn0 plt.plot(bestCN0, bestAcc, 'bo') - # 1/12T, Tcoh=20 ms - elif mode == "doperr" and doperr <= 1.0 / (12 * 0.02): + # 1/12T, Tcoh=20 ms, Tcoh=5 ms + elif mode == "doperr" and doperr <= 1.0 / (12 * Tcoh): bestAcc = acc bestCN0 = cn0 plt.plot(bestCN0, bestAcc, 'bo') @@ -122,12 +145,15 @@ def dynamicAccPlot(filename, mode): #plt.plot(cn0VecAll, accVecAll, 'bo') plt.plot(bestCN0, bestAcc, 'r.-') plt.xlabel('CN0') - plt.ylabel('Acceleration Hz/s') - plt.grid() + plt.ylabel('Acceleration (g)') + plt.grid(True) if mode == "lockrate": plt.title("PLL lock rate >= 0.68 (1-sigma)") elif mode == "doperr": - plt.title("Doppler 1-sigma error <= 1/(12*0.02) = 4.2 Hz") + if cfg.isL1CA(): + plt.title("Doppler 1-sigma error <= 1/(12*0.005) = 16.7 Hz") + else: + plt.title("Doppler 1-sigma error <= 1/(12*0.02) = 4.2 Hz") plt.show() @@ -147,8 +173,12 @@ def main(): parser.add_argument("-e", "--dyn-acc-dop", help="x-sigma Doppler error acceleration tolerance vs. CN0", action="store_true") + parser.add_argument("-b", "--band", + help="l1ca or l2c (default)") args = parser.parse_args() + if args.band: + cfg.BAND = args.band if args.lockrate: lockrateCn0Plot(args.filename) elif args.dyn_acc_lockrate: diff --git a/peregrine/analysis/run_sim.py b/peregrine/analysis/run_sim.py index 9929cb7..a5a6dc1 100755 --- a/peregrine/analysis/run_sim.py +++ b/peregrine/analysis/run_sim.py @@ -20,81 +20,153 @@ import peregrine.iqgen.bits.signals as signals import numpy as np -IQ_DATA = "iqdata.bin" -TRACK_DATA = "track_res" -TRACK_RES_DATA = "track_res.json" + +class CfgClass: + + def __init__(self): + self.IQ_DATA = "iqdata.bin" + self.TRACK_DATA = "track_res" + self.TRACK_RES_DATA = "track_res.json" + self.BAND = "l2c" # "l1ca" + self.fpgaSim = " " # "--short-long-cycles " + + def isL1CA(self): + if "l1ca" == self.BAND: + return True + return False + + def __str__(self): + s = "band:" + self.BAND + "\n" + s = s + "fpga delay control simulation: " + self.fpgaSim + return s + + +cfg = CfgClass() def runCmd(cmd): print cmd # Do not use shell=True for untrusted input! - out = subprocess.check_output(cmd, shell=True) + try: + out = subprocess.check_output(cmd, shell=True) + except subprocess.CalledProcessError as e: + # Peregrine exits with failure if acquisition fails. + # Therefore handle error case + print e.output + return e.output, False + print out - return out + return out, True def runIqGen(lens, snr, dop, acc): cmd = "python " + peregrinePath() + "/peregrine/iqgen/iqgen_main.py" - cmd = cmd + " --gps-sv 1 --encoder 1bit --bands l1ca+l2c --message-type crc --profile low_rate" + if cfg.isL1CA(): + cmd = cmd + " --gps-sv 1 --encoder 1bit --bands l1ca --message-type crc --profile low_rate" + else: + cmd = cmd + " --gps-sv 1 --encoder 1bit --bands l1ca+l2c --message-type crc --profile low_rate" + if float(acc) != 0.0: cmd = cmd + " --doppler-type linear --doppler-speed " + acc + " " elif float(dop) != 0.0: cmd = cmd + " --doppler-type const " + dop + " " else: cmd = cmd + " --doppler-type zero " - cmd = cmd + " --snr " + snr + " --generate " + lens + " " - cmd = cmd + " --output " + IQ_DATA + " " - out = runCmd(cmd) + cmd = cmd + "--amplitude-type poly --amplitude-units snr-db --amplitude-a0 " + snr + cmd = cmd + " --generate " + lens + " " + + cmd = cmd + " --output " + cfg.IQ_DATA + " " + out, success = runCmd(cmd) + if not success: + print "iqgen failed" + sys.exit(1) + lines = out.split('\n') cn0 = None - l2dop = None - l2chip = None + l_dop = None + l_chip = None for ln in lines: words = ln.split() if len(words) == 0: continue - if words[0] == ".L2": - cn0 = words[2] - if words[0] == ".l2_doppler:": - l2dop = words[1] - if words[0] == ".l2_chip:": - l2chip = words[1] + if cfg.isL1CA(): + if words[0] == ".l1": + cn0 = words[2] + if words[0] == ".l1_doppler:": + l_dop = words[1] + if words[0] == ".l1_chip:": + l_chip = words[1] + else: + if words[0] == ".l2": + cn0 = words[2] + if words[0] == ".l2_doppler:": + l_dop = words[1] + if words[0] == ".l2_chip:": + l_chip = words[1] + if words[0] == ".SNR": if float(words[2]) != float(snr): print "snr unexpected" sys.exit(1) - if cn0 is None or l2dop is None or l2chip is None: + if cn0 is None or l_dop is None or l_chip is None: print "iqgen output parse error" sys.exit(1) - return lens, snr, l2dop, acc, l2chip, cn0, cmd + return lens, snr, l_dop, acc, l_chip, cn0, cmd def runTracker(dopp, cp, lens): - cmd = "python " + peregrinePath() + \ - "/peregrine/analysis/tracking_loop.py -f 1bit_x2 -P 1 --profile low_rate " - cmd = cmd + "-p " + cp + " -d " + dopp - cmd = cmd + " -o " + TRACK_DATA + " " - cmd = cmd + " -S l2c " - cmd = cmd + IQ_DATA + " " - out = runCmd(cmd) + if cfg.isL1CA(): + cmd = "python " + peregrinePath() + \ + "/peregrine/analysis/tracking_loop.py -f 1bit -P 1 --profile low_rate --l1ca-profile med --ms-to-process -1 " + # "--short-long-cycles " tracking loop corrections are taken in use in FPGA with a delay of 1 ms + cmd = cmd + cfg.fpgaSim + cmd = cmd + "-p " + cp + " -d " + dopp + cmd = cmd + " -o " + cfg.TRACK_DATA + " " + cmd = cmd + " -S l1ca " + cmd = cmd + " --file " + cfg.IQ_DATA + " " + + else: + cmd = "python " + peregrinePath() + \ + "/peregrine/analysis/tracking_loop.py -f 1bit_x2 -P 1 --profile low_rate --ms-to-process -1 " + cmd = cmd + "-p " + cp + " -d " + dopp + cmd = cmd + " -o " + cfg.TRACK_DATA + " " + cmd = cmd + " -S l2c " + cmd = cmd + " --file " + cfg.IQ_DATA + " " + + out, success = runCmd(cmd) lines = out.split('\n') + # This is usefull if acquisition is run instead of + # running tracking_loop.py directly + if not success: + for ln in lines: + if ln.find("No satellites acquired"): + # Acceptable failure + return cmd, False + # Unacceptable failure + print "Acquisition/tracking failed unexpectedly" + sys.exit(1) + durationOk = False for ln in lines: words = ln.split() if len(words) == 0: continue if words[0] == "Time" and words[1] == "to" and words[2] == "process": - if int(words[4]) / 1000 == int(lens): + if round(float(words[4])) == int(lens): durationOk = True if not durationOk: print "Data duration mismatch" sys.exit(1) - return cmd + return cmd, True def processTrackResults(acc): - data = np.genfromtxt(TRACK_DATA + ".PRN-1.l2c", - dtype=float, delimiter=',', names=True) + if cfg.isL1CA(): + data = np.genfromtxt(cfg.TRACK_DATA + ".PRN-1.l1ca", + dtype=float, delimiter=',', names=True) + else: + data = np.genfromtxt(cfg.TRACK_DATA + ".PRN-1.l2c", + dtype=float, delimiter=',', names=True) CN0 = data['CN0'] dopp = data['carr_doppler'] lock = data['lock_detect_outp'] @@ -105,10 +177,15 @@ def processTrackResults(acc): lockRate = np.sum(lock) / nsamples dopErr = np.ndarray(shape=(1, len(dopp)), dtype=float) - i = np.linspace(1, len(dopp), len(dopp), dtype=float) - # Doppler - (i * 20 ms * acc Hz/s * L2/L1 Hz relation) - dopErr = np.abs(dopp - i * 0.02 * acc * - (signals.GPS.L2C.CENTER_FREQUENCY_HZ / signals.GPS.L1CA.CENTER_FREQUENCY_HZ)) + + if cfg.isL1CA(): + ms_tracked = data['ms_tracked'] + dopErr = np.abs(dopp - ms_tracked / 1000.0 * acc) + else: + i = np.linspace(1, len(dopp), len(dopp), dtype=float) + # Doppler - (i * 20 ms * acc Hz/s * L2/L1 Hz relation) + dopErr = np.abs(dopp - i * 0.02 * acc * + (signals.GPS.L2C.CENTER_FREQUENCY_HZ / signals.GPS.L1CA.CENTER_FREQUENCY_HZ)) maxDopErr = np.max(dopErr) sortedDopErr = np.sort(dopErr) ix1 = int(0.5 + 0.6827 * (len(sortedDopErr) - 1)) @@ -126,33 +203,35 @@ def produce(lens, snr, dop, acc): snr = str(float(snr)) dop = str(float(dop)) acc = str(float(acc)) - lens, snr, l2dop, acc, l2chip, cn0, iqgenCmd = runIqGen(lens, snr, dop, acc) - trackerCmd = runTracker(l2dop, l2chip, lens) - avgCN0, lockRate, dopSigma1, dopSigma2, dopSigma3, maxDopErr = processTrackResults( - acc) - fpout = open(TRACK_RES_DATA, "a") - d = datetime.datetime(2000, 1, 1) - js = json.dumps({"stamp": d.utcnow().isoformat(), - "snr": snr, - "duration": float(lens), - "iqgencn0": float(cn0), - "l2dop": float(l2dop), - "l2chip": float(l2chip), - "acc": float(acc), - "avgcn0": avgCN0, - "lockrate": lockRate, - "dopSigma1": dopSigma1, - "dopSigma2": dopSigma2, - "dopSigma3": dopSigma3, - "dopMaxErr": maxDopErr, - "iqgencmd": iqgenCmd, - "trackcmd": trackerCmd}, - sort_keys=True, indent=4, separators=(',', ': ')) - print js - fpout.write(js + ',') - fpout.close() - - return lockRate, dopSigma1, dopSigma2 + lens, snr, l_dop, acc, l_chip, cn0, iqgenCmd = runIqGen(lens, snr, dop, acc) + trackerCmd, trackSuccess = runTracker(l_dop, l_chip, lens) + if trackSuccess: + avgCN0, lockRate, dopSigma1, dopSigma2, dopSigma3, maxDopErr = processTrackResults( + acc) + fpout = open(cfg.TRACK_RES_DATA, "a") + d = datetime.datetime(2000, 1, 1) + js = json.dumps({"stamp": d.utcnow().isoformat(), + "snr": snr, + "duration": float(lens), + "iqgencn0": float(cn0), + "l_dop": float(l_dop), + "l_chip": float(l_chip), + "acc": float(acc), + "avgcn0": avgCN0, + "lockrate": lockRate, + "dopSigma1": dopSigma1, + "dopSigma2": dopSigma2, + "dopSigma3": dopSigma3, + "dopMaxErr": maxDopErr, + "iqgencmd": iqgenCmd, + "trackcmd": trackerCmd}, + sort_keys=True, indent=4, separators=(',', ': ')) + print js + fpout.write(js + ',') + fpout.close() + return lockRate, dopSigma1, dopSigma2, trackSuccess + else: + return 0, 0, 0, trackSuccess def runCn0Range(): @@ -161,7 +240,7 @@ def runCn0Range(): doppler = 0 # Hz acceleration = 0.0 # Hz / s for snr in snrRng: - lockRate, dopSigma1, dopSigma2 = produce( + lockRate, dopSigma1, dopSigma2, success = produce( length, snr / 10.0, doppler, acceleration) @@ -188,10 +267,10 @@ def runDynamicLockRate(): bestLockRate = lockRate bestAcc = acc - lockRate, dopSigma1, dopSigma2 = produce( + lockRate, dopSigma1, dopSigma2, success = produce( length, snr / 10.0, doppler, acc) print "SNR", snr / 10.0, "ACC", acc, "LOCKRATE", lockRate - if lockRate >= lockRateThreshold: + if lockRate >= lockRateThreshold and success: acc = acc + accStep else: print "BEST ACC", bestAcc, "LOCKRATE", bestLockRate @@ -205,7 +284,12 @@ def runDynamicFreq(): # Must run from low CN0 to high CN0 doppler = 0 # Hz accStep = 10 # Size of single acceleration step (L1CA Hz/s) - freqErrThreshold = 1.0 / (12 * 0.02) # 1/12T, Tcoh=20 ms. Failure threshold. + if cfg.isL1CA(): + # 1/12T, Tcoh=5 ms. Failure threshold. + freqErrThreshold = 1.0 / (12 * 0.005) + else: + # 1/12T, Tcoh=20 ms. Failure threshold. + freqErrThreshold = 1.0 / (12 * 0.02) # When it is reached, the next CN0 bin is tried freqErr = 0.0 # 1-sigma @@ -220,11 +304,11 @@ def runDynamicFreq(): bestFreqErr = freqErr bestAcc = acc - lockRate, dopSigma1, dopSigma2 = produce( + lockRate, dopSigma1, dopSigma2, trackSuccess = produce( length, snr / 10.0, doppler, acc) freqErr = dopSigma1 print "SNR", snr / 10.0, "ACC", acc, "1-SIGMA", freqErr - if freqErr <= freqErrThreshold: + if freqErr <= freqErrThreshold and trackSuccess: acc = acc + accStep else: print "BEST ACC", bestAcc, "1-SIGMA", bestFreqErr @@ -237,23 +321,32 @@ def peregrinePath(): def main(): - global TRACK_RES_DATA parser = argparse.ArgumentParser() parser.add_argument("-l", "--lockrate", help="Simple lockrate vs. CN0", action="store_true") parser.add_argument("-f", "--filename", - help="Output file which is appended. Default " + TRACK_RES_DATA) + help="Output file which is appended. Default " + cfg.TRACK_RES_DATA) parser.add_argument("-d", "--dyn-lockrate", help="Lockrate, acceleration, CN0", action="store_true") parser.add_argument("-e", "--dyn-freq", help="Fequency error, acceleration, CN0", action="store_true") + parser.add_argument("-b", "--band", + help="l1ca or l2c (default)") + parser.add_argument("-s", "--short-long-cycles", + help="FPGA delay control simulation", + action="store_true") args = parser.parse_args() if args.filename: - TRACK_RES_DATA = args.filename + cfg.TRACK_RES_DATA = args.filename + if args.band: + cfg.BAND = args.band + if args.short_long_cycles: + cfg.fpgaSim = "--short-long-cycles " + if args.lockrate: runCn0Range() elif args.dyn_lockrate: