CMakeified almost everything. Test code in python/ and apps other than uhd_modes.py still need minor updating.

apps/CMakeLists.txt

@@ -0,0 +1,26 @@
+# Copyright 2011 Free Software Foundation, Inc.
+# 
+# This file is part of GNU Radio
+# 
+# GNU Radio is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 3, or (at your option)
+# any later version.
+# 
+# GNU Radio is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+# 
+# You should have received a copy of the GNU General Public License
+# along with GNU Radio; see the file COPYING.  If not, write to
+# the Free Software Foundation, Inc., 51 Franklin Street,
+# Boston, MA 02110-1301, USA.
+
+include(GrPython)
+
+GR_PYTHON_INSTALL(
+    PROGRAMS
+    uhd_modes.py
+    DESTINATION bin
+)

apps/get_correlated_records.py

@@ -0,0 +1,132 @@
+#!/usr/bin/env python
+from modes_parse import modes_parse
+import mlat
+import numpy
+import sys
+
+#sffile = open("27augsf3.txt")
+#rudifile = open("27augrudi3.txt")
+
+#sfoutfile = open("sfout.txt", "w")
+#rudioutfile = open("rudiout.txt", "w")
+
+sfparse = modes_parse([37.762236,-122.442525])
+
+sf_station = [37.762236,-122.442525, 100]
+mv_station = [37.409348,-122.07732, 100]
+bk_station = [37.854246, -122.266701, 100]
+
+raw_stamps = []
+
+#first iterate through both files to find the estimated time difference. doesn't have to be accurate to more than 1ms or so.
+#to do this, look for type 17 position packets with the same data. assume they're unique. print the tdiff.
+
+#collect a list of raw timestamps for each aircraft from each station
+#the raw stamps have to be processed into corrected stamps OR distance has to be included in each
+#then postprocess to find clock delay for each and determine drift rate for each aircraft separately
+#then come up with an average clock drift rate
+#then find an average drift-corrected clock delay
+#then find rms error
+
+#ok so get [ICAO, [raw stamps], [distance]] for each matched record
+
+files = [open(arg) for arg in sys.argv[1:]]
+
+#files = [sffile, rudifile]
+stations = [sf_station, mv_station]#, bk_station]
+
+records = []
+
+for each_file in files:
+    recordlist = []
+    for line in each_file:
+        [msgtype, shortdata, longdata, parity, ecc, reference, timestamp] = line.split()
+        recordlist.append({"data": {"msgtype": long(msgtype, 10),\
+                                    "shortdata": long(shortdata, 16),\
+                                    "longdata": long(longdata, 16),\
+                                    "parity": long(parity, 16),\
+                                    "ecc": long(ecc, 16)},
+                           "time": float(timestamp)\
+                          })
+    records.append(recordlist)
+
+#ok now we have records parsed into something usable that we can == with
+
+def feet_to_meters(feet):
+    return feet * 0.3048006096012
+
+all_heard = []
+#gather list of reports which were heard by all stations
+for station0_report in records[0]: #iterate over list of reports from station 0
+    for other_reports in records[1:]:
+        stamps = [station0_report["time"]]
+        stamp = [report["time"] for report in other_reports if report["data"] == station0_report["data"]]# for other_reports in records[1:]]
+        if len(stamp) > 0:
+            stamps.append(stamp[0])
+    if len(stamps) == len(records): #found same report in all records
+        all_heard.append({"data": station0_report["data"], "times": stamps})
+
+#print all_heard
+
+#ok, now let's pull out the location-bearing packets so we can find our time offset
+position_reports = [x for x in all_heard if x["data"]["msgtype"] == 17 and 9 <= (x["data"]["longdata"] >> 51) & 0x1F <= 18]
+offset_list = []
+#there's probably a way to list-comprehension-ify this but it looks hard
+for msg in position_reports:
+    data = msg["data"]
+    [alt, lat, lon, rng, bearing] = sfparse.parseBDS05(data["shortdata"], data["longdata"], data["parity"], data["ecc"])
+    ac_pos = [lat, lon, feet_to_meters(alt)]
+    rel_times = []
+    for time, station in zip(msg["times"], stations):
+        #here we get the estimated time at the aircraft when it transmitted
+        range_to_ac = numpy.linalg.norm(numpy.array(mlat.llh2ecef(station))-numpy.array(mlat.llh2ecef(ac_pos)))
+        timestamp_at_ac = time - range_to_ac / mlat.c
+        rel_times.append(timestamp_at_ac)
+    offset_list.append({"aircraft": data["shortdata"] & 0xffffff, "times": rel_times})
+
+#this is a list of unique aircraft, heard by all stations, which transmitted position packets
+#we do drift calcs separately for each aircraft in the set because mixing them seems to screw things up
+#i haven't really sat down and figured out why that is yet
+unique_aircraft = list(set([x["aircraft"] for x in offset_list]))
+print "Aircraft heard for clock drift estimate: %s" % [str("%x" % ac) for ac in unique_aircraft]
+print "Total reports used: %d over %.2f seconds" % (len(position_reports), position_reports[-1]["times"][0]-position_reports[0]["times"][0])
+
+#get a list of reported times gathered by the unique aircraft that transmitted them
+#abs_unique_times = [report["times"] for ac in unique_aircraft for report in offset_list if report["aircraft"] == ac]
+#print abs_unique_times
+#todo: the below can probably be done cleaner with nested list comprehensions
+clock_rate_corrections = [0]
+for i in range(1,len(stations)):
+    drift_error_limited = []
+    for ac in unique_aircraft:
+        times = [report["times"] for report in offset_list if report["aircraft"] == ac]
+
+        s0_times = [report[0] for report in times]
+        rel_times = [report[i]-report[0] for report in times]
+
+        #find drift error rate
+        drift_error = [(y-x)/(b-a) for x,y,a,b in zip(rel_times, rel_times[1:], s0_times[0:], s0_times[1:])]
+        drift_error_limited.append([x for x in drift_error if abs(x) < 1e-5])
+
+    #flatten the list of lists (tacky, there's a better way)
+    drift_error_limited = [x for sublist in drift_error_limited for x in sublist]
+    clock_rate_corrections.append(0-numpy.mean(drift_error_limited))
+
+for i in range(len(clock_rate_corrections)):
+    print "drift from %d relative to station 0: %.3fppm" % (i, clock_rate_corrections[i] * 1e6)
+
+#let's get the average clock offset (based on drift-corrected, TDOA-corrected derived timestamps)
+clock_offsets = [[numpy.mean([x["times"][i]*(1+clock_rate_corrections[i])-x["times"][0] for x in offset_list])][0] for i in range(0,len(stations))]
+for i in range(len(clock_offsets)):
+    print "mean offset from %d relative to station 0: %.3f seconds" % (i, clock_offsets[i])
+
+#for the two-station case, let's now go back, armed with our clock drift and offset, and get the variance between expected and observed timestamps
+error_list = []
+for i in range(1,len(stations)):
+    for report in offset_list:
+        error = abs(((report["times"][i]*(1+clock_rate_corrections[i]) - report["times"][0]) - clock_offsets[i]) * mlat.c)
+        error_list.append(error)
+        #print error
+
+rms_error = (numpy.mean([error**2 for error in error_list]))**0.5
+print "RMS error in TDOA: %.1f meters" % rms_error

apps/mlat-2station-postanalysis.py

@@ -0,0 +1,40 @@
+#!/usr/bin/env python
+import numpy
+import mlat
+
+#rudi says:
+#17 8da12615 903bf4bd3eb2c0 36ac95 000000 0.0007421782357 2.54791875
+#17 8d4b190a 682de4acf8c177 5b8f55 000000 0.0005142348236 2.81227225
+
+#sf says:
+#17 8da12615 903bf4bd3eb2c0 36ac95 000000 0.003357535461 00.1817445
+#17 8d4b190a 682de4acf8c177 5b8f55 000000 0.002822938375 000.446215
+
+sf_station = [37.762236,-122.442525, 100]
+mv_station = [37.409348,-122.07732, 100]
+
+report1_location = [37.737804, -122.485139, 3345]
+report1_sf_tstamp = 0.1817445
+report1_mv_tstamp = 2.54791875
+
+report2_location = [37.640836, -122.260218, 2484]
+report2_sf_tstamp = 0.446215
+report2_mv_tstamp = 2.81227225
+
+report1_tof_sf = numpy.linalg.norm(numpy.array(mlat.llh2ecef(sf_station))-numpy.array(mlat.llh2ecef(report1_location))) / mlat.c
+report1_tof_mv = numpy.linalg.norm(numpy.array(mlat.llh2ecef(mv_station))-numpy.array(mlat.llh2ecef(report1_location))) / mlat.c
+
+report1_sf_tstamp_abs = report1_sf_tstamp - report1_tof_sf
+report1_mv_tstamp_abs = report1_mv_tstamp - report1_tof_mv
+
+report2_tof_sf = numpy.linalg.norm(numpy.array(mlat.llh2ecef(sf_station))-numpy.array(mlat.llh2ecef(report2_location))) / mlat.c
+report2_tof_mv = numpy.linalg.norm(numpy.array(mlat.llh2ecef(mv_station))-numpy.array(mlat.llh2ecef(report2_location))) / mlat.c
+
+report2_sf_tstamp_abs = report2_sf_tstamp - report2_tof_sf
+report2_mv_tstamp_abs = report2_mv_tstamp - report2_tof_mv
+
+dt1 = report1_sf_tstamp_abs - report1_mv_tstamp_abs
+dt2 = report2_sf_tstamp_abs - report2_mv_tstamp_abs
+
+error = abs((dt1-dt2) * mlat.c)
+print error

apps/uhd_modes.py

@@ -0,0 +1,204 @@
+#!/usr/bin/env python
+# Copyright 2010 Nick Foster
+# 
+# This file is part of gr-air-modes
+# 
+# gr-air-modes is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 3, or (at your option)
+# any later version.
+# 
+# gr-air-modes is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+# 
+# You should have received a copy of the GNU General Public License
+# along with gr-air-modes; see the file COPYING.  If not, write to
+# the Free Software Foundation, Inc., 51 Franklin Street,
+# Boston, MA 02110-1301, USA.
+# 
+
+my_position = [37.76225, -122.44254]
+#my_position = [37.409066,-122.077836]
+#my_position = None
+
+from gnuradio import gr, gru, optfir, eng_notation, blks2
+from gnuradio import uhd
+from gnuradio.eng_option import eng_option
+from optparse import OptionParser
+import time, os, sys, threading
+from string import split, join
+import air_modes
+import gnuradio.gr.gr_threading as _threading
+import csv
+
+class top_block_runner(_threading.Thread):
+    def __init__(self, tb):
+        _threading.Thread.__init__(self)
+        self.setDaemon(1)
+        self.tb = tb
+        self.done = False
+        self.start()
+
+    def run(self):
+        self.tb.run()
+        self.done = True
+
+class adsb_rx_block (gr.top_block):
+  def __init__(self, options, args, queue):
+    gr.top_block.__init__(self)
+
+    self.options = options
+    self.args = args
+    rate = int(options.rate)
+
+    if options.filename is None:
+      self.u = uhd.single_usrp_source("", uhd.io_type_t.COMPLEX_FLOAT32, 1)
+      time_spec = uhd.time_spec(0.0)
+      self.u.set_time_now(time_spec)
+
+      #if(options.rx_subdev_spec is None):
+      #  options.rx_subdev_spec = ""
+      #self.u.set_subdev_spec(options.rx_subdev_spec)
+      if not options.antenna is None:
+        self.u.set_antenna(options.antenna)
+
+      self.u.set_samp_rate(rate)
+      rate = int(self.u.get_samp_rate()) #retrieve actual
+
+      if options.gain is None: #set to halfway
+        g = self.u.get_gain_range()
+        options.gain = (g.start()+g.stop()) / 2.0
+
+      if not(self.tune(options.freq)):
+        print "Failed to set initial frequency"
+
+      print "Setting gain to %i" % (options.gain,)
+      self.u.set_gain(options.gain)
+      print "Gain is %i" % (self.u.get_gain(),)
+
+    else:
+      self.u = gr.file_source(gr.sizeof_gr_complex, options.filename)
+
+    print "Rate is %i" % (rate,)
+
+    pass_all = 0
+    if options.output_all :
+      pass_all = 1
+
+    self.demod = gr.complex_to_mag()
+    self.avg = gr.moving_average_ff(100, 1.0/100, 400)
+    
+    #the DBSRX especially tends to be spur-prone; the LPF keeps out the
+    #spur multiple that shows up at 2MHz
+    self.lpfiltcoeffs = gr.firdes.low_pass(1, rate, 1.8e6, 100e3)
+    self.lpfilter = gr.fir_filter_ccf(1, self.lpfiltcoeffs)
+    
+    self.preamble = air_modes.modes_preamble(rate, options.threshold)
+    #self.framer = air_modes.modes_framer(rate)
+    self.slicer = air_modes.modes_slicer(rate, queue)
+    
+    self.connect(self.u, self.lpfilter, self.demod)
+    self.connect(self.demod, self.avg)
+    self.connect(self.demod, (self.preamble, 0))
+    self.connect(self.avg, (self.preamble, 1))
+    self.connect((self.preamble, 0), (self.slicer, 0))
+
+  def tune(self, freq):
+    result = self.u.set_center_freq(freq, 0)
+    return result
+
+def printraw(msg):
+    print msg
+
+if __name__ == '__main__':
+  usage = "%prog: [options] output filename"
+  parser = OptionParser(option_class=eng_option, usage=usage)
+  parser.add_option("-R", "--rx-subdev-spec", type="string",
+          help="select USRP Rx side A or B", metavar="SUBDEV")
+  parser.add_option("-A", "--antenna", type="string",
+          help="select which antenna to use on daughterboard")
+  parser.add_option("-f", "--freq", type="eng_float", default=1090e6,
+                      help="set receive frequency in Hz [default=%default]", metavar="FREQ")
+  parser.add_option("-g", "--gain", type="int", default=None,
+                      help="set RF gain", metavar="dB")
+  parser.add_option("-r", "--rate", type="eng_float", default=4000000,
+                      help="set ADC sample rate [default=%default]")
+  parser.add_option("-T", "--threshold", type="eng_float", default=3.0,
+                      help="set pulse detection threshold above noise in dB [default=%default]")
+  parser.add_option("-a","--output-all", action="store_true", default=False,
+                      help="output all frames")
+  parser.add_option("-F","--filename", type="string", default=None,
+            help="read data from file instead of USRP")
+  parser.add_option("-K","--kml", type="string", default=None,
+                      help="filename for Google Earth KML output")
+  parser.add_option("-P","--sbs1", action="store_true", default=False,
+                      help="open an SBS-1-compatible server on port 30003")
+  parser.add_option("-w","--raw", action="store_true", default=False,
+                      help="open a server outputting raw timestamped data on port 9988")
+  parser.add_option("-n","--no-print", action="store_true", default=False,
+                      help="disable printing decoded packets to stdout")
+  parser.add_option("-l","--location", type="string", default=None,
+                      help="GPS coordinates of receiving station in format xx.xxxxx,xx.xxxxx")
+  (options, args) = parser.parse_args()
+
+  if options.location is not None:
+    reader = csv.reader([options.location], quoting=csv.QUOTE_NONNUMERIC)
+    my_position = reader.next()
+
+  queue = gr.msg_queue()
+  
+  outputs = [] #registry of plugin output functions
+  updates = [] #registry of plugin update functions
+
+  if options.kml is not None:
+    sqlport = air_modes.modes_output_sql(my_position, 'adsb.db') #create a SQL parser to push stuff into SQLite
+    outputs.append(sqlport.insert)
+    #also we spawn a thread to run every 30 seconds (or whatever) to generate KML
+    kmlgen = modes_kml('adsb.db', options.kml, my_position) #create a KML generating thread which reads the database
+
+  if options.sbs1 is True:
+    sbs1port = air_modes.modes_output_sbs1(my_position)
+    outputs.append(sbs1port.output)
+    updates.append(sbs1port.add_pending_conns)
+    
+  if options.no_print is not True:
+    outputs.append(air_modes.modes_output_print(my_position).parse)
+
+  if options.raw is True:
+    rawport = air_modes.modes_raw_server()
+    outputs.append(rawport.output)
+    outputs.append(printraw)
+    updates.append(rawport.add_pending_conns)
+
+  fg = adsb_rx_block(options, args, queue)
+  runner = top_block_runner(fg)
+
+  while 1:
+    try:
+      #the update registry is really for the SBS1 and raw server plugins -- we're looking for new TCP connections.
+      #i think we have to do this here rather than in the output handler because otherwise connections will stack up
+      #until the next output arrives
+      for update in updates:
+        update()
+      
+      #main message handler
+      if queue.empty_p() == 0 :
+        while queue.empty_p() == 0 :
+          msg = queue.delete_head() #blocking read
+
+          for out in outputs:
+            out(msg.to_string())
+
+      elif runner.done:
+        raise KeyboardInterrupt
+      else:
+        time.sleep(0.1)
+
+    except KeyboardInterrupt:
+      fg.stop()
+      runner = None
+      if options.kml is not None:
+          kmlgen.done = True
+      break