from collections import namedtuple, defaultdict import sqlite3 from datetime import datetime import time import rpyc import numpy as np from scipy.misc import pilutil import scipy.ndimage.filters as image_filters import scipy.ndimage as ndimage from utils.path import path from utils import getPyPlot fmtts = lambda ts: str(datetime.fromtimestamp(ts)).split('.')[0] UTC_OFFSET = 3600 * 2 # 2 hours summertime WISPY_OFFSET = 0.7 CHANNELS = range(1, 15) BURST_PERIOD = 6.0 ZERO_PERIOD = 6.0 CHANNELS_B = { #~ 0: (2401, 2495), 1: (2401, 2423), 2: (2406, 2428), 3: (2411, 2433), 4: (2416, 2438), 5: (2421, 2443), 6: (2426, 2448), 7: (2431, 2453), 8: (2436, 2458), 9: (2441, 2463), 10: (2446, 2468), 11: (2451, 2473), 12: (2456, 2478), 13: (2461, 2483), 14: (2473, 2495) } #~ CHANNELS_B = dict([(e[0], (e[1][0] + 7, e[1][1] - 7)) for e in CHANNELS_B.items()]) BURST_STATIONS = { '107': ('127.0.0.1', 18107), '108': ('127.0.0.1', 18108), '110': ('127.0.0.1', 18110), } def runBurst(station): ''' connect to AP and send bursts over different channels''' host, port = BURST_STATIONS[station] channel2ts = {} conn = rpyc.classic.connect(host, port=port) conn.modules.sys.path.append('/sbin') chaburst = conn.modules.chaburst try: chaburst.stopDefaultWifi() chaburst.startBurstDevice() for chan in CHANNELS: print 'measure channel %s' % chan chaburst.setChannel(chan) ts1 = time.time() chaburst.sendBurst(BURST_PERIOD) channel2ts[chan] = (ts1, time.time()) finally: chaburst.startDefaultWifi() conn.close() print 'measure zero channel' # store timestamps for a period without traffic at fake channel 0 #~ ts1 = time.time() #~ time.sleep(ZERO_PERIOD) #~ channel2ts[0] = (ts1, time.time()) return channel2ts def gatherTimestampsForPositions(stations): ''' collect measurements at different locations interactivly ''' print 'Enter first location index:', loc2ts = {} # key (station, location), value {channel: (ts_start, ts_end)} while True: try: loc = raw_input() loc = int(loc) except ValueError: if loc == 'q': print 'done...' break else: print 'need a int value for location!' continue ts1 = time.time() for station in stations: print 'start measure %s for location %s at station %s' % (fmtts(ts1), loc, station) #~ print 'end measurement with [return]', #~ raw_input() try: channel2ts = runBurst(station) except Exception, KeyboardInterrupt: print 'error during running burst' continue #~ ts2 = time.time() loc2ts[(station, loc)] = channel2ts print 'continue with next location: ', print 'got the following locations' for (loc, station), ts in loc2ts.items(): print loc, station, ts return loc2ts def getSignalStrength(data): '''estimate the paramters of a 1d GMM with 2 modes''' from scikits.learn.mixture import GMM # assume 2 modes - one gaussian for the noise, and one for the signal clf = GMM(n_states=2, cvtype='full') clf.fit(np.array(data).reshape(len(data), 1)) # assume the larger mean is our signal signalStrengthMax = max(clf.means) idx = list(clf.means).index(signalStrengthMax) varianceMax = list(clf.covars)[idx] signalStrengthMin = min(clf.means) return signalStrengthMax, varianceMax, signalStrengthMin def updateSignalData(locations, ssid, wispy_db, loc2ts): ''' update Location() data in locations with measured rssi values for given ssid the values are extracted from the wispy_db file by using timestamps from measurements_file ''' # currently only first ts is used for locid, (ts1, ts2) in loc2ts.items(): loc2ts[locid] = ts1 # now use these timestamps to extract signals from wispy-db firstts, lastts = min(loc2ts.values()), max(loc2ts.values()) res = {} with sqlite3.connect(wispy_db) as conn: c = conn.cursor() sql = '''select rssi, milliseconds_since_epoch from network_data where milliseconds_since_epoch > ? and milliseconds_since_epoch < ? and ssid = ?''' c.execute(sql, ((firstts - 10 + UTC_OFFSET) * 1000, (lastts + 10 + UTC_OFFSET) * 1000, ssid)) ts2rssi = {} for row in c: ts = row[1] / 1000 - UTC_OFFSET ts2rssi[ts] = row[0] #~ print datetime.fromtimestamp(ts) for loc, ts in loc2ts.items(): # iterate over all stored rssi values # and pick that one with the minimum timestamp delta (thats not fast but should be ok) delttats, key = min([(abs(ts - e), e) for e in ts2rssi.keys()]) #~ print '%s at location %s with rssi: %s [delta %.2f sec]' % (fmtts(ts), loc, ts2rssi[key], delttats) locations[loc] = locations[loc]._replace(rssi=ts2rssi[key], ts=ts) def ts2dbts(ts): ''' convert ts to milliseconds_since_epoch ''' return (ts + UTC_OFFSET) * 1000 def dbts2ts(dbts): return dbts / 1000.0 - UTC_OFFSET def rasterize(xs, ys, zs, readings_per_sweep, destfile): first_ts = np.min(xs) last_ts = np.max(xs) delta = last_ts - first_ts #~ print delta resolution_x = len(set(xs)) resolution_y = readings_per_sweep data = np.zeros((resolution_x, resolution_y)) for x, y, z in zip(xs, ys, zs): _x = int(round((x - first_ts) / float(delta) * (resolution_x - 1))) data[_x][y] = z data = pilutil.imresize(data, (1000, resolution_y)) data = data.transpose() pilutil.imsave(destfile, data) Location = namedtuple('Location', 'id x y z rssi ts') def loadLocations(f): print 'loading locations from %s' % f.abspath() lines = [l for l in f.lines() if l.strip() and not l.strip().startswith('#')] res = {} for l in lines: i, x, y, z = [float(e.strip()) for e in l.split()] i = int(i) res[i] = Location(i, x, y, z, None, None) return res def loadMeasurements(): ''' return a list of locations with measured rssi signal strengths ''' p = path(r'D:\loco-dev\maps\whs19\experiment\measure_locs.txt') locs = loadLocations(p) wispy_db = path('d:/loco-dev/maps/whs19/experiment/measure.wsx') run_file = path('d:/loco-dev/maps/whs19/experiment/run1.txt') updateSignalData(locs, 'CS31', wispy_db, run_file) return locs Sweep = namedtuple('Sweep', 'dt data setting channel') Setting = namedtuple('Setting', 'id starting_frequency frequency_resolution readings_per_sweep amplitude_offset amplitude_resolution rssi_max') class Measure(object): def __init__(self, location_file, wispy_db, plotpath, loc2tsfile): self.location_file = location_file self.wispy_db = wispy_db self.loc2tsfile = loc2tsfile if not self.loc2tsfile.exists(): self.loc2tsfile.touch() # contents wispy_db setting table self.settings = {} # key: setting id, value: Setting() # load data into self.settings self.loadWiSpySettings() self.loc2ts = defaultdict(list) # key: location id, values: list of dicts containing timestamps for each channel # load data into self.loc2ts self.loadMeasurements() self.plotpath = plotpath if not plotpath.exists(): plotpath.makedirs() def loadWiSpySettings(self): self.settings.clear() with sqlite3.connect(self.wispy_db) as conn: c = conn.cursor() # load all settings sql = '''select id, starting_frequency, frequency_resolution, readings_per_sweep, amplitude_offset, amplitude_resolution, rssi_max from setting ''' c.execute(sql) for row in c: self.settings[row[0]] = Setting(*row) def loadSweepDataForLocation(self, locid, station): sweeps = defaultdict(list) # keys: measure_idx, values: [Sweep(), ] with sqlite3.connect(self.wispy_db) as conn: c = conn.cursor() channeldata = self.loc2ts[(locid, station)] for measure_idx, channel2ts in enumerate(channeldata): allts = [] for ts1, ts2 in channel2ts.values(): allts.extend([ts1, ts2]) fromts, untilts = min(allts), max(allts) # load sweep data from firsts to lastts sql = '''select sweep_data, milliseconds_since_epoch, setting_id from sweep where milliseconds_since_epoch > ? and milliseconds_since_epoch < ? order by milliseconds_since_epoch ''' MORE_SECONDS = 10 # fetch more seconds at the left/right side of intervall c.execute(sql, ((fromts - MORE_SECONDS + UTC_OFFSET) * 1000, (untilts + MORE_SECONDS + UTC_OFFSET) * 1000)) for row in c: ts = dbts2ts(row[1]) for channel, (ts1, ts2) in channel2ts.items(): if ts1 <= ts <= ts2: sweep = Sweep(datetime.fromtimestamp(ts), row[0], row[2], channel) sweeps[measure_idx].append(sweep) # check if we found any data in wispy db if not measure_idx in sweeps: print datetime.fromtimestamp(fromts - MORE_SECONDS) print 'got no sweep data from wispy db for measurement %s at location %s' % (measure_idx, locid) #~ print sweeps return sweeps def combineWispyDBs(self, source, dest, clear=False): if clear: with sqlite3.connect(dest) as conn: c = conn.cursor() c.execute('delete from network_data') c.execute('delete from sweep') c.execute('delete from setting') with sqlite3.connect(source) as conn: c = conn.cursor() # select last setting sql = '''select id, device, start_time_msepoch, starting_frequency, frequency_resolution, readings_per_sweep, amplitude_offset, amplitude_resolution, rssi_max from setting order by id desc''' setting = list(c.execute(sql))[-1] sql = '''select milliseconds_since_epoch, sweep_data from sweep where setting_id = ?''' sweeps = list(c.execute(sql, str(setting[0]))) print 'got %s sweeps' % len(sweeps) sql = '''select ssid, channel, encryption, mac, rssi, supported_rates, signal_quality, mode, milliseconds_since_epoch from network_data''' network_data = list(c.execute(sql)) print 'inserting data into %s' % dest with sqlite3.connect(dest) as conn: c = conn.cursor() sql = '''select max(id) from setting''' try: setting_id = list(c.execute(sql))[0][0] + 1 except Exception: setting_id = 0 print 'using new setting id %s' % setting_id sql = '''insert into setting (id, device, start_time_msepoch, starting_frequency, frequency_resolution, readings_per_sweep, amplitude_offset, amplitude_resolution, rssi_max) values(?, ?, ?, ?, ?, ?, ?, ?, ?)''' c.execute(sql, [setting_id] + list(setting[1:])) sql = '''select max(id) from sweep''' try: next_sweep_id = list(c.execute(sql))[0][0] + 1 except Exception: next_sweep_id = 0 sql = '''insert into sweep (id, setting_id, milliseconds_since_epoch, sweep_data) values(?, ?, ?, ?)''' for s_id, sweep in enumerate(sweeps, next_sweep_id): c.execute(sql, [s_id, setting_id] + list(sweep)) sql = '''select max(id) from network_data''' try: next_nd_id = list(c.execute(sql))[0][0] + 1 except Exception: next_nd_id = 0 sql = '''insert into network_data (id, ssid, channel, encryption, mac, rssi, supported_rates, signal_quality, mode, milliseconds_since_epoch) values(?, ?, ?, ?, ?, ? , ?, ?, ?, ?)''' for nd_id, nd in enumerate(network_data, next_nd_id): c.execute(sql, [nd_id] + list(nd)) def runMeasurement(self, stations): ## check if stations are reachable for s in stations: try: conn = rpyc.classic.connect(BURST_STATIONS[s][0], port=BURST_STATIONS[s][1]) except: raise RuntimeError('station %s is not reachable' % s) else: print 'station %s is reachable' % s conn.close() loc2ts = gatherTimestampsForPositions(stations) s = '\n'.join('%s/%s %r' % (station, locid, channel2ts) for (station, locid), channel2ts in loc2ts.items()) print 'storing measurements result to %s' % self.loc2tsfile.abspath() self.loc2tsfile.write_text('run from %s at locations %s\n%s\n\n' % (datetime.now(), ','.join(str(e) for e in loc2ts), s), append=True) # refresh self.loc2ts self.loadMeasurements() self.loadWiSpySettings() print 'please save wispy db now!', raw_input() def loadMeasurements(self): ''' load measurements from loc2tsfile (next to location_file) and return dict with keys: location id, values: list of dicts that contain (start/end) timestamps for each channel ''' self.loc2ts.clear() print 'loading measurements from %s' % self.loc2tsfile.abspath() for line in [l for l in self.loc2tsfile.lines() if l.strip()]: if line.startswith('run from'): continue #line format: 107/1 {1: (1315234649.809, 1315234651.031), 2: (1315234651 ... l = line.split() station, locid = l[0].split('/') locid = int(locid) channel2ts = eval(' '.join(l[1:])) self.loc2ts[(locid, station)].append(channel2ts) # adjust wispy timestamp offset (empirically determined!) for channeldata in self.loc2ts.values(): for channel2ts in channeldata: for chanid, (ts1, ts2) in channel2ts.items(): channel2ts[chanid] = (ts1 + WISPY_OFFSET, ts2 + WISPY_OFFSET) def analyzeMeasurements(self, locid, measure_idx, station): #~ locs = loadLocations(self.location_file) sweep = self.loadSweepDataForLocation(locid, station) self.drawMeasure(sweep[measure_idx], locid, station, measure_idx) avg_sig = self.analyzeBurst(sweep[measure_idx], locid, measure_idx, draw=True) print 'got avg sig %.2f' % avg_sig def buildOptimizeFile(self, optimizefile, station): measure_idx = 0 locs = loadLocations(self.location_file) loc2sig = [] for locid, loc in locs.items(): sweep = self.loadSweepDataForLocation(locid, station) avg_sig = self.analyzeBurst(sweep[measure_idx], locid, measure_idx) if avg_sig is not None: loc2sig.append((loc, avg_sig)) s = '\n'.join('%s %s %s %s %.2f' % (loc.id, loc.x, loc.y, loc.z, avg_sig) for loc, avg_sig in loc2sig) optimizefile = path(optimizefile % station) print 'storing optimize file to %s' % optimizefile.abspath() optimizefile.write_text(s) def drawMeasure(self, sweeps, locid, station, measure_idx): channel2ts = self.loc2ts[(locid, station)][measure_idx] plt, mlab, dates, font_manager, ticker = getPyPlot() width = 1200 aspectRatio = 0.5 dpi = 70 figsize = width / float(dpi), width / float(dpi) * aspectRatio fig = plt.figure(figsize=figsize) ax1 = plt.axes([0.08, 0.21, 0.85, 0.75]) # assert all Sweep() instances are from the same wispy setting assert len(set([s.setting for s in sweeps])) == 1 setting = self.settings[sweeps[0].setting] starting_frequency = setting.starting_frequency frequency_resolution = setting.frequency_resolution / 1000 amplitude_offset = setting.amplitude_offset amplitude_resolution = setting.amplitude_resolution readings_per_sweep = setting.readings_per_sweep xs, ys, zs = [], [], [] for sweep in sweeps: assert len(sweep.data) == readings_per_sweep for i, c in enumerate(sweep.data): # TODO: is this correct? amplitude = amplitude_offset + ord(c) * amplitude_resolution #~ print amplitude if amplitude > -100: xs.append(sweep.dt) y = (starting_frequency + (i * frequency_resolution)) / 1000.0 #~ print starting_frequency, i, frequency_resolution ys.append(i) zs.append(ord(c)) xs = mlab.date2num(xs) _ys = np.array(ys) _ys = ((_ys * frequency_resolution) + starting_frequency) / 1000 ax1.scatter(xs, _ys , s=2, c='#AAAACC', linewidths=0) for channel, (ts1, ts2) in channel2ts.items(): _xs = mlab.date2num([datetime.fromtimestamp(ts) for ts in (ts1, ts2)]) ax1.plot(_xs, [CHANNELS_B[channel]] * 2, c='#7766FF') ax1.xaxis_date() ax1.xaxis.set_major_formatter(dates.DateFormatter('%H:%M:%S')) fig.savefig(self.plotpath.joinpath('full_%s_%s.png' % (locid, measure_idx)), format='png', dpi=dpi) plt.close(fig) rasterize(xs, ys, zs, readings_per_sweep, self.plotpath.joinpath('rast_%s_%s.png' % (locid, measure_idx))) def default_analyzer(self, channeldata): sigstrenghts = {} for c, cdata in channeldata.items(): if c == 0: continue sigstrenghts[c] = getSignalStrength(cdata) MIN_STRENGTH = 3 # ignore signals that are too near at sig0 all_sigs = [sig for sig, variance, sig0 in sigstrenghts.values() if sig > sig0 + MIN_STRENGTH and variance < 50] if len(all_sigs) > 0: avg_sig = sum(all_sigs) / float(len(all_sigs)) else: avg_sig = -100 return sigstrenghts, avg_sig def analyzeBurst(self, sweeps, locid, measure_idx, draw=False): c = len(set([s.setting for s in sweeps])) if c != 1: print 'anticipated exactly one setting for locid %s, measure_idx %s, but got %s' % (locid, measure_idx, c) return setting = self.settings[sweeps[0].setting] #~ y = (starting_frequency + (i * frequency_resolution)) / 1000.0 channeldata = defaultdict(list) for sweep in sweeps: minfreq , maxfreq = CHANNELS_B[sweep.channel] _toidx = lambda f: int((f - (setting.starting_frequency / 1000.0)) / (setting.frequency_resolution / 10.0**6)) minidx = _toidx(minfreq) maxidx = _toidx(maxfreq) for c in sweep.data[minidx:maxidx]: ampl = setting.amplitude_offset + ord(c) * setting.amplitude_resolution channeldata[sweep.channel].append(ampl) sigstrenghts, avg_sig = default_analyzer(channeldata) if draw: plt, mlab, dates, font_manager, ticker = getPyPlot() for c, cdata in channeldata.items(): if c == 0: continue fig = plt.figure(figsize=(12, 6)) ax = plt.axes() r, BIN_WIDTH = (-120, 0), 2 bins = [e for e in range(r[0], r[1], BIN_WIDTH)] ax.hist(cdata, bins=bins) ax.set_xlim(r) sig, variance, sig0 = sigstrenghts[c] ax.axvline(sig) ax.axvline(sig0, color='red') ax.axvline(avg_sig, color='green', linewidth=2.0) ax1 = ax.twinx() ax1.plot(range(-120, 50), [mlab.normpdf(x, sig, variance)[0][0] for x in range(-120, 50)]) ax1.set_xlim(r) plt.title('loc %s, run %s, chan %s [mu: %.1f, sigma: %.1f, mu0: %.1f, mu_total: %.2f]' % (locid, measure_idx, c, sig, variance, sig0, avg_sig)) fig.savefig(self.plotpath.joinpath('channel_%02d_%02d_%02d.png' % (locid, measure_idx, c)), dpi=80) plt.close(fig) return avg_sig if __name__ == '__main__': #python /usr/lib/python2.6/site-packages/rpyc/scripts/rpyc_classic.py -m forking -p 20001 #~ wispy_db = path('d:/loco-dev/maps/whs19/experiment/measure.wsx') #~ location_file = path(r'D:\loco-dev\maps\whs19\experiment\measure_locs.txt') location_file = path(r'D:\loco-dev\maps\umic\experiment\locations.txt') wispy_db = path('d:/loco-dev/maps/umic/experiment/measure.wsx') loc2tsfile = path(r'D:\loco-dev\maps\umic\experiment\loc2ts.txt') plotpath = path('r:/') optimizefile = path('d:/loco-dev/maps/umic/experiment/measure_%s.txt') m = Measure(location_file, wispy_db, plotpath, loc2tsfile) stations = ['107', '108', '110'] #~ stations = ['110'] #~ stations = ['110'] #~ m.runMeasurement(stations) m.combineWispyDBs('d:/loco-dev/maps/umic/experiment/measure1.wsx', 'd:/loco-dev/maps/umic/experiment/measure.wsx', clear=True) m.combineWispyDBs('d:/loco-dev/maps/umic/experiment/measure2.wsx', 'd:/loco-dev/maps/umic/experiment/measure.wsx') m.combineWispyDBs('d:/loco-dev/maps/umic/experiment/measure3.wsx', 'd:/loco-dev/maps/umic/experiment/measure.wsx') m.combineWispyDBs('d:/loco-dev/maps/umic/experiment/measure5.wsx', 'd:/loco-dev/maps/umic/experiment/measure.wsx',) #~ m.combineWispyDBs('d:/loco-dev/maps/umic/experiment/measure4.wsx', 'd:/loco-dev/maps/umic/experiment/measure.wsx') #~ m.combineWispyDBs('d:/loco-dev/maps/umic/experiment/measure6.wsx', 'd:/loco-dev/maps/umic/experiment/measure.wsx') #~ for i in range(1, 12): #~ m.analyzeMeasurements(i, 0) #~ m.analyzeMeasurements(3, 0, '110') #~ m.analyzeMeasurements(4, 0, '110') #~ m.analyzeMeasurements(15, 0, '110') for s in stations: m.buildOptimizeFile(optimizefile, s)