From 3c8c60f57fb691028f868c159713db90bb219482 Mon Sep 17 00:00:00 2001 From: Nick Foster Date: Sun, 14 Oct 2012 17:56:01 -0700 Subject: [PATCH 1/4] Removing emitter-centered local decoding from the CPR parser. Using only global decoding. --- python/cpr.py | 27 +++++---------------------- python/msprint.py | 4 +++- 2 files changed, 8 insertions(+), 23 deletions(-) diff --git a/python/cpr.py b/python/cpr.py index 5115bf5..1d647dd 100755 --- a/python/cpr.py +++ b/python/cpr.py @@ -172,7 +172,6 @@ def range_bearing(loc_a, loc_b): class cpr_decoder: def __init__(self, my_location): self.my_location = my_location - self.lkplist = {} self.evenlist = {} self.oddlist = {} @@ -180,9 +179,9 @@ class cpr_decoder: self.my_location = new_location def weed_poslists(self): - for poslist in [self.lkplist, self.evenlist, self.oddlist]: + for poslist in [self.evenlist, self.oddlist]: for key, item in poslist.items(): - if time.time() - item[2] > 900: + if time.time() - item[2] > 10: del poslist[key] def decode(self, icao24, encoded_lat, encoded_lon, cpr_format, surface): @@ -194,30 +193,14 @@ class cpr_decoder: [decoded_lat, decoded_lon] = [None, None] - #okay, let's traverse the lists and weed out those entries that are older than 15 minutes, as they're unlikely to be useful. + #okay, let's traverse the lists and weed out those entries that are older than 10 seconds self.weed_poslists() - if icao24 in self.lkplist: - #do emitter-centered local decoding - [decoded_lat, decoded_lon] = cpr_resolve_local(self.lkplist[icao24][0:2], [encoded_lat, encoded_lon], cpr_format, surface) - self.lkplist[icao24] = [decoded_lat, decoded_lon, time.time()] #update the local position for next time - - elif (icao24 in self.evenlist) \ - and (icao24 in self.oddlist) \ - and (abs(self.evenlist[icao24][2] - self.oddlist[icao24][2]) < 10) \ - and (surface == 0): + if (icao24 in self.evenlist) \ + and (icao24 in self.oddlist): newer = (self.oddlist[icao24][2] - self.evenlist[icao24][2]) > 0 #figure out which report is newer [decoded_lat, decoded_lon] = cpr_resolve_global(self.evenlist[icao24][0:2], self.oddlist[icao24][0:2], newer, surface) #do a global decode - self.lkplist[icao24] = [decoded_lat, decoded_lon, time.time()] - #so we really can't guarantee that local decoding will work unless you are POSITIVE that you can't hear more than 180nm out. - #this will USUALLY work, but you can't guarantee it! - #elif surface == 1 and self.my_location is not None: - # [local_lat, local_lon] = cpr_resolve_local(self.my_location, [encoded_lat, encoded_lon], cpr_format, surface) #try local decoding - # [rnge, bearing] = range_bearing(self.my_location, [local_lat, local_lon]) - # if rnge < validrange: #if the local decoding can be guaranteed valid - # self.lkplist[icao24] = [local_lat, local_lon, time.time()] #update the local position for next time - # [decoded_lat, decoded_lon] = [local_lat, local_lon] else: raise CPRNoPositionError diff --git a/python/msprint.py b/python/msprint.py index d613fa4..2349c35 100644 --- a/python/msprint.py +++ b/python/msprint.py @@ -69,7 +69,9 @@ class output_print(air_modes.parse): pass def output(self, msg): - print self.parse(msg) + parsed = self.parse(msg) + if parsed is not None: + print self.parse(msg) def print0(self, shortdata, ecc): [vs, cc, sl, ri, altitude] = self.parse0(shortdata) From ac9e09893b118c90e4f1fca66cd9287cce3fd01f Mon Sep 17 00:00:00 2001 From: Nick Foster Date: Sun, 14 Oct 2012 22:34:50 -0700 Subject: [PATCH 2/4] CPR changes. Fixed a lot of the surface position errors. At least one bug still remains. --- python/cpr.py | 89 ++++++++++++++++++++++++++++++--------------------- 1 file changed, 52 insertions(+), 37 deletions(-) diff --git a/python/cpr.py b/python/cpr.py index 1d647dd..7c29013 100755 --- a/python/cpr.py +++ b/python/cpr.py @@ -28,12 +28,6 @@ from air_modes.exceptions import * latz = 15 -def nbits(surface): - if surface == 1: - return 19 - else: - return 17 - def nz(ctype): return 4 * latz - ctype @@ -49,9 +43,6 @@ def dlat(ctype, surface): else: return tmp / nzcalc -def nl_eo(declat_in, ctype): - return nl(declat_in) - ctype - def nl(declat_in): if abs(declat_in) >= 87.0: return 1.0 @@ -62,25 +53,23 @@ def dlon(declat_in, ctype, surface): tmp = 90.0 else: tmp = 360.0 - nlcalc = nl_eo(declat_in, ctype) + nlcalc = nl(declat_in)-ctype if nlcalc == 0: return tmp else: - return tmp / max(nlcalc, 1.0) + return tmp / nlcalc def decode_lat(enclat, ctype, my_lat, surface): tmp1 = dlat(ctype, surface) - tmp2 = float(enclat) / (2**nbits(surface)) + tmp2 = float(enclat) / (2**17) j = math.floor(my_lat/tmp1) + math.floor(0.5 + ((my_lat % tmp1) / tmp1) - tmp2) -# print "dlat gives " + "%.6f " % tmp1 + "with j = " + "%.6f " % j + " and tmp2 = " + "%.6f" % tmp2 + " given enclat " + "%x" % enclat return tmp1 * (j + tmp2) def decode_lon(declat, enclon, ctype, my_lon, surface): tmp1 = dlon(declat, ctype, surface) - tmp2 = float(enclon) / (2.0**nbits(surface)) + tmp2 = float(enclon) / (2**17) m = math.floor(my_lon / tmp1) + math.floor(0.5 + ((my_lon % tmp1) / tmp1) - tmp2) -# print "dlon gives " + "%.6f " % tmp1 + "with m = " + "%.6f " % m + " and tmp2 = " + "%.6f" % tmp2 + " given enclon " + "%x" % enclon return tmp1 * (m + tmp2) @@ -93,7 +82,11 @@ def cpr_resolve_local(my_location, encoded_location, ctype, surface): return [decoded_lat, decoded_lon] -def cpr_resolve_global(evenpos, oddpos, mostrecent, surface): +def cpr_resolve_global(evenpos, oddpos, mypos, mostrecent, surface): + #cannot resolve surface positions unambiguously without knowing receiver position + if surface and mypos is None: + raise CPRNoPositionError + dlateven = dlat(0, surface) dlatodd = dlat(1, surface) @@ -114,7 +107,6 @@ def cpr_resolve_global(evenpos, oddpos, mostrecent, surface): #This checks to see if the latitudes of the reports straddle a transition boundary #If so, you can't get a globally-resolvable location. if nl(rlateven) != nl(rlatodd): - #print "Boundary straddle!" raise CPRBoundaryStraddleError if mostrecent == 0: @@ -122,21 +114,41 @@ def cpr_resolve_global(evenpos, oddpos, mostrecent, surface): else: rlat = rlatodd + #disambiguate latitude + if surface: + if mypos[0] < 0: + rlat -= 90 + dl = dlon(rlat, mostrecent, surface) nlthing = nl(rlat) - ni = max(nlthing - mostrecent, 1) + ni = nlthing - mostrecent - m = math.floor(((evenpos[1]*(nlthing-1)-oddpos[1]*(nlthing))/2**17)+0.5) #longitude index + #THIS COLLAPSES WHEN M IS NOT A FACTOR OF FOUR + m = math.floor(((evenpos[1]*(nlthing-1)-oddpos[1]*nlthing)/2**17)+0.5) #longitude index + if m%4: + print "IM GONNA DIE" + + print "DL: %f nlthing: %f ni: %f m: %f" % (dl, nlthing, ni, m) if mostrecent == 0: enclon = evenpos[1] else: enclon = oddpos[1] - rlon = dl * (((ni+m) % ni)+enclon/2**17) + rlon = dl * ((m % ni)+enclon/2**17) + if surface: + #longitudes need to be resolved to the nearest 90 degree segment to the receiver. + wat = mypos[1] + if wat < 0: + wat += 360 + print "mylon: %f mylon unwrapped: %f rlon raw: %f" % (mypos[1], wat, rlon) + zone = lambda lon: 90 * (int(lon) / 90) + rlon += (zone(wat) - zone(rlon)) + + #limit to (-180, 180) if rlon > 180: - rlon = rlon - 360.0 + rlon -= 360.0 return [rlat, rlon] @@ -199,8 +211,7 @@ class cpr_decoder: if (icao24 in self.evenlist) \ and (icao24 in self.oddlist): newer = (self.oddlist[icao24][2] - self.evenlist[icao24][2]) > 0 #figure out which report is newer - [decoded_lat, decoded_lon] = cpr_resolve_global(self.evenlist[icao24][0:2], self.oddlist[icao24][0:2], newer, surface) #do a global decode - + [decoded_lat, decoded_lon] = cpr_resolve_global(self.evenlist[icao24][0:2], self.oddlist[icao24][0:2], self.my_location, newer, surface) #do a global decode else: raise CPRNoPositionError @@ -223,11 +234,11 @@ def cpr_encode(lat, lon, ctype, surface): yz = math.floor(scalar * ((lat % dlati)/dlati) + 0.5) rlat = dlati * ((yz / scalar) + math.floor(lat / dlati)) - nleo = nl_eo(rlat, ctype) + nleo = nl(rlat)-ctype if nleo == 0: dloni = 360.0 else: - dloni = 360.0 / nl_eo(rlat, ctype) + dloni = 360.0 / nleo xz = math.floor(scalar * ((lon % dloni)/dloni) + 0.5) @@ -245,53 +256,57 @@ if __name__ == '__main__': #latitudes the corresponding error in longitude is greater bs = 0 + surface = True for i in range(0, rounds): - decoder = cpr_decoder(None) even_lat = random.uniform(-85, 85) - even_lon = random.uniform(-180,180) + even_lon = random.uniform(-180, 180) odd_lat = even_lat + 1e-2 odd_lon = min(even_lon + 1e-2, 180) + decoder = cpr_decoder([odd_lat, odd_lon]) #encode that position - (evenenclat, evenenclon) = cpr_encode(even_lat, even_lon, False, False) - (oddenclat, oddenclon) = cpr_encode(odd_lat, odd_lon, True, False) + (evenenclat, evenenclon) = cpr_encode(even_lat, even_lon, False, surface) + (oddenclat, oddenclon) = cpr_encode(odd_lat, odd_lon, True, surface) #perform a global decode icao = random.randint(0, 0xffffff) try: - evenpos = decoder.decode(icao, evenenclat, evenenclon, False, False) + evenpos = decoder.decode(icao, evenenclat, evenenclon, False, surface) #print "CPR global decode with only one report: %f %f" % (evenpos[0], evenpos[1]) raise Exception("CPR test failure: global decode with only one report") except CPRNoPositionError: pass try: - (odddeclat, odddeclon, rng, brg) = decoder.decode(icao, oddenclat, oddenclon, True, False) + (odddeclat, odddeclon, rng, brg) = decoder.decode(icao, oddenclat, oddenclon, True, surface) except CPRBoundaryStraddleError: bs += 1 continue except CPRNoPositionError: raise Exception("CPR test failure: no decode after even/odd inputs") - #print "Lat: %f Lon: %f" % (ac_lat, ac_lon) - if abs(odddeclat - odd_lat) > threshold or abs(odddeclon - odd_lon) > threshold: - print "odddeclat: %f odd_lat: %f" % (odddeclat, odd_lat) - print "odddeclon: %f odd_lon: %f" % (odddeclon, odd_lon) + print "F odddeclat: %f odd_lat: %f" % (odddeclat, odd_lat) + print "F odddeclon: %f odd_lon: %f" % (odddeclon, odd_lon) raise Exception("CPR test failure: global decode error greater than threshold") + else: + print "S odddeclat: %f odd_lat: %f" % (odddeclat, odd_lat) + print "S odddeclon: %f odd_lon: %f" % (odddeclon, odd_lon) nexteven_lat = odd_lat + 1e-2 nexteven_lon = min(odd_lon + 1e-2, 180) - (nexteven_enclat, nexteven_enclon) = cpr_encode(nexteven_lat, nexteven_lon, False, False) + (nexteven_enclat, nexteven_enclon) = cpr_encode(nexteven_lat, nexteven_lon, False, surface) try: - (evendeclat, evendeclon) = cpr_resolve_local([even_lat, even_lon], [nexteven_enclat, nexteven_enclon], False, False) + (evendeclat, evendeclon) = cpr_resolve_local([even_lat, even_lon], [nexteven_enclat, nexteven_enclon], False, surface) except CPRNoPositionError: raise Exception("CPR test failure: local decode failure to resolve") if abs(evendeclat - nexteven_lat) > threshold or abs(evendeclon - nexteven_lon) > threshold: + print "F evendeclat: %f nexteven_lat: %f evenlat: %f" % (evendeclat, nexteven_lat, even_lat) + print "F evendeclon: %f nexteven_lon: %f evenlon: %f" % (evendeclon, nexteven_lon, even_lon) raise Exception("CPR test failure: local decode error greater than threshold") print "CPR test successful. There were %i boundary straddles over %i rounds." % (bs, rounds) From 0b30db353daab873eb8add9588a237495b97263a Mon Sep 17 00:00:00 2001 From: Nick Foster Date: Wed, 17 Oct 2012 09:13:10 -0700 Subject: [PATCH 3/4] CPR: surface decoding still exhibits occasional errors when m~=0. Suspect it might actually be a bug in the encoder. --- python/cpr.py | 93 +++++++++++++++++++++++++++------------------------ 1 file changed, 50 insertions(+), 43 deletions(-) diff --git a/python/cpr.py b/python/cpr.py index 7c29013..decefc0 100755 --- a/python/cpr.py +++ b/python/cpr.py @@ -26,6 +26,8 @@ from air_modes.exceptions import * #the decoder is implemented as a class, cpr_decoder, which keeps state for local decoding. #the encoder is cpr_encode([lat, lon], type (even=0, odd=1), and surface (0 for surface, 1 for airborne)) +#TODO: remove range/bearing calc from CPR decoder class. you can do this outside of the decoder. + latz = 15 def nz(ctype): @@ -46,18 +48,15 @@ def dlat(ctype, surface): def nl(declat_in): if abs(declat_in) >= 87.0: return 1.0 - return math.floor( (2.0*math.pi) * pow(math.acos(1.0- (1.0-math.cos(math.pi/(2.0*latz))) / pow( math.cos( (math.pi/180.0)*abs(declat_in) ) ,2.0) ),-1.0)) + return math.floor( (2.0*math.pi) * math.acos(1.0- (1.0-math.cos(math.pi/(2.0*latz))) / math.cos( (math.pi/180.0)*abs(declat_in) )**2 )**-1) def dlon(declat_in, ctype, surface): - if surface == 1: + if surface: tmp = 90.0 else: tmp = 360.0 - nlcalc = nl(declat_in)-ctype - if nlcalc == 0: - return tmp - else: - return tmp / nlcalc + nlcalc = max(nl(declat_in)-ctype, 1) + return tmp / nlcalc def decode_lat(enclat, ctype, my_lat, surface): tmp1 = dlat(ctype, surface) @@ -120,29 +119,29 @@ def cpr_resolve_global(evenpos, oddpos, mypos, mostrecent, surface): rlat -= 90 dl = dlon(rlat, mostrecent, surface) - nlthing = nl(rlat) - ni = nlthing - mostrecent + nl_rlat = nl(rlat) - #THIS COLLAPSES WHEN M IS NOT A FACTOR OF FOUR - m = math.floor(((evenpos[1]*(nlthing-1)-oddpos[1]*nlthing)/2**17)+0.5) #longitude index - if m%4: - print "IM GONNA DIE" - - print "DL: %f nlthing: %f ni: %f m: %f" % (dl, nlthing, ni, m) + m = math.floor(((evenpos[1]*(nl_rlat-1)-oddpos[1]*nl_rlat)/2**17)+0.5) #longitude index + + #when surface positions straddle a disambiguation boundary (90 degrees), + #surface decoding will fail. this might never be a problem in real life, but it'll fail in the + #test case. the documentation doesn't mention it. if mostrecent == 0: enclon = evenpos[1] else: enclon = oddpos[1] - rlon = dl * ((m % ni)+enclon/2**17) + rlon = dl * ((m % max(nl_rlat-mostrecent,1)) + enclon/2.**17) + + #print "DL: %f nl: %f m: %f rlon: %f" % (dl, nl_rlat, m, rlon) + #print "evenpos: %x, oddpos: %x, mostrecent: %i" % (evenpos[1], oddpos[1], mostrecent) if surface: #longitudes need to be resolved to the nearest 90 degree segment to the receiver. wat = mypos[1] if wat < 0: wat += 360 - print "mylon: %f mylon unwrapped: %f rlon raw: %f" % (mypos[1], wat, rlon) zone = lambda lon: 90 * (int(lon) / 90) rlon += (zone(wat) - zone(rlon)) @@ -226,58 +225,64 @@ class cpr_decoder: #encode CPR position def cpr_encode(lat, lon, ctype, surface): if surface is True: - scalar = float(2**19) + scalar = 2.**19 else: - scalar = float(2**17) + scalar = 2.**17 - dlati = float(dlat(ctype, False)) + #encode using 360 constant for segment size. + dlati = dlat(ctype, False) yz = math.floor(scalar * ((lat % dlati)/dlati) + 0.5) rlat = dlati * ((yz / scalar) + math.floor(lat / dlati)) - nleo = nl(rlat)-ctype - if nleo == 0: - dloni = 360.0 - else: - dloni = 360.0 / nleo - + #encode using 360 constant for segment size. + dloni = dlon(lat, ctype, False) xz = math.floor(scalar * ((lon % dloni)/dloni) + 0.5) - yz = int(yz % scalar) - xz = int(xz % scalar) + yz = int(yz) & (2**17-1) + xz = int(xz) & (2**17-1) return (yz, xz) #lat, lon if __name__ == '__main__': import sys, random - - rounds = 10000 + + rounds = 10001 threshold = 1e-3 #0.001 deg lat/lon #this accuracy is highly dependent on latitude, since at high #latitudes the corresponding error in longitude is greater bs = 0 - surface = True + surface = False + + lats = [i/(rounds/170.)-85 for i in range(0,rounds)] + lons = [i/(rounds/360.)-180 for i in range(0,rounds)] for i in range(0, rounds): - even_lat = random.uniform(-85, 85) - even_lon = random.uniform(-180, 180) - odd_lat = even_lat + 1e-2 - odd_lon = min(even_lon + 1e-2, 180) + even_lat = lats[i] + #even_lat = random.uniform(-85, 85) + even_lon = lons[i] + #even_lon = random.uniform(-180, 180) + odd_lat = even_lat + 1e-3 + odd_lon = min(even_lon + 1e-3, 180) decoder = cpr_decoder([odd_lat, odd_lon]) #encode that position (evenenclat, evenenclon) = cpr_encode(even_lat, even_lon, False, surface) (oddenclat, oddenclon) = cpr_encode(odd_lat, odd_lon, True, surface) - #perform a global decode + #try to perform a global decode -- this should fail since the decoder + #only has heard one position. need two for global decoding. icao = random.randint(0, 0xffffff) try: evenpos = decoder.decode(icao, evenenclat, evenenclon, False, surface) - #print "CPR global decode with only one report: %f %f" % (evenpos[0], evenpos[1]) raise Exception("CPR test failure: global decode with only one report") except CPRNoPositionError: pass + #now try to do a real decode with the last packet's odd complement + #watch for a boundary straddle -- this isn't fatal, it just indicates + #that the even and odd reports lie on either side of a longitudinal boundary + #and so you can't get a position try: (odddeclat, odddeclon, rng, brg) = decoder.decode(icao, oddenclat, oddenclon, True, surface) except CPRBoundaryStraddleError: @@ -290,20 +295,22 @@ if __name__ == '__main__': print "F odddeclat: %f odd_lat: %f" % (odddeclat, odd_lat) print "F odddeclon: %f odd_lon: %f" % (odddeclon, odd_lon) raise Exception("CPR test failure: global decode error greater than threshold") - else: - print "S odddeclat: %f odd_lat: %f" % (odddeclat, odd_lat) - print "S odddeclon: %f odd_lon: %f" % (odddeclon, odd_lon) +# else: +# print "S odddeclat: %f odd_lat: %f" % (odddeclat, odd_lat) +# print "S odddeclon: %f odd_lon: %f" % (odddeclon, odd_lon) - nexteven_lat = odd_lat + 1e-2 - nexteven_lon = min(odd_lon + 1e-2, 180) + nexteven_lat = odd_lat + 1e-3 + nexteven_lon = min(odd_lon + 1e-3, 180) (nexteven_enclat, nexteven_enclon) = cpr_encode(nexteven_lat, nexteven_lon, False, surface) + #try a locally-referenced decode try: (evendeclat, evendeclon) = cpr_resolve_local([even_lat, even_lon], [nexteven_enclat, nexteven_enclon], False, surface) except CPRNoPositionError: raise Exception("CPR test failure: local decode failure to resolve") - + + #check to see if the positions were valid if abs(evendeclat - nexteven_lat) > threshold or abs(evendeclon - nexteven_lon) > threshold: print "F evendeclat: %f nexteven_lat: %f evenlat: %f" % (evendeclat, nexteven_lat, even_lat) print "F evendeclon: %f nexteven_lon: %f evenlon: %f" % (evendeclon, nexteven_lon, even_lon) From a3f2dec1f935856c25af58a8c6ef8b177749fb73 Mon Sep 17 00:00:00 2001 From: Nick Foster Date: Wed, 17 Oct 2012 09:24:15 -0700 Subject: [PATCH 4/4] CPR decoder: Keep a separate list for sfc and air positions to keep them from stepping on each other. --- python/cpr.py | 25 +++++++++++++++++++------ 1 file changed, 19 insertions(+), 6 deletions(-) diff --git a/python/cpr.py b/python/cpr.py index decefc0..c16289b 100755 --- a/python/cpr.py +++ b/python/cpr.py @@ -185,6 +185,8 @@ class cpr_decoder: self.my_location = my_location self.evenlist = {} self.oddlist = {} + self.evenlist_sfc = {} + self.oddlist_sfc = {} def set_location(new_location): self.my_location = new_location @@ -194,23 +196,34 @@ class cpr_decoder: for key, item in poslist.items(): if time.time() - item[2] > 10: del poslist[key] + for poslist in [self.evenlist_sfc, self.oddlist_sfc]: + for key, item in poslist.items(): + if time.time() - item[2] > 25: + del poslist[key] def decode(self, icao24, encoded_lat, encoded_lon, cpr_format, surface): + if surface: + oddlist = self.oddlist_sfc + evenlist = self.evenlist_sfc + else: + oddlist = self.oddlist + evenlist = self.evenlist + #add the info to the position reports list for global decoding if cpr_format==1: - self.oddlist[icao24] = [encoded_lat, encoded_lon, time.time()] + oddlist[icao24] = [encoded_lat, encoded_lon, time.time()] else: - self.evenlist[icao24] = [encoded_lat, encoded_lon, time.time()] + evenlist[icao24] = [encoded_lat, encoded_lon, time.time()] [decoded_lat, decoded_lon] = [None, None] #okay, let's traverse the lists and weed out those entries that are older than 10 seconds self.weed_poslists() - if (icao24 in self.evenlist) \ - and (icao24 in self.oddlist): - newer = (self.oddlist[icao24][2] - self.evenlist[icao24][2]) > 0 #figure out which report is newer - [decoded_lat, decoded_lon] = cpr_resolve_global(self.evenlist[icao24][0:2], self.oddlist[icao24][0:2], self.my_location, newer, surface) #do a global decode + if (icao24 in evenlist) \ + and (icao24 in oddlist): + newer = (oddlist[icao24][2] - evenlist[icao24][2]) > 0 #figure out which report is newer + [decoded_lat, decoded_lon] = cpr_resolve_global(evenlist[icao24][0:2], oddlist[icao24][0:2], self.my_location, newer, surface) #do a global decode else: raise CPRNoPositionError