PK*J44clocklogger/__init__.py"""Clock and weather logger""" __version__ = '0.1' PK[*J$v,,clocklogger/analysis.pyfrom __future__ import division import numpy as np from numpy import sin, cos, pi from datetime import timedelta PLOT = False try: import matplotlib.pyplot as plt except ImportError: PLOT = False if PLOT: plt.ioff() PPS_RATE_ERROR_THRESHOLD = 50e-6 class DataError(Exception): pass def round_time_to_three_seconds(t): return t.replace(second = (t.second//3)*3, microsecond = 0) class ClockAnalyser(object): def __init__(self, source, initial_drift=0): self.source = source self.drift_offset = initial_drift self.last_drift = None self.cache = np.empty((0, 2)) self.cache_start_time = None self.pretrigger = 0.2 # seconds self.edge_level = 0.2 self.shim_width = 24.0 # millimetres ~= milliradians (at 1m from pivot) self.decay_fit_duration = 0.05 # duration of fit segment (s) self.decay_fit_delay = 0.003 # wait after crossing threshold (s) self.decay_fit_level = 0.5 # y-value to use as reference point def process(self, pps_edge='up', sampling_rate_from_pps=False, fit_decay=False): """Read samples from source and yield drift & amplitude values""" #pretrigger_samples = self.pretrigger * source.fs last_time = None for t, (i_pos_pps, i_neg_pps), (i_pos_tick, i_neg_tick) in self.generate_edge_groups(fit_decay=fit_decay): i_pps = i_pos_pps if pps_edge == 'up' else i_neg_pps self.sanity_check_pps(i_pps) # Find first PPS after down tick & calculate drift drift_delta = self.calculate_drift(i_pos_tick, i_pps, sampling_rate_from_pps) # Calculate amplitude amplitude = self.calculate_amplitude(i_pos_tick, i_neg_tick) print(("Found ticks: PPS %5d up, %5d down\n" + " Tick %5d up, %5d down") \ % tuple([len(_) for _ in [i_pos_pps, i_neg_pps, i_pos_tick, i_neg_tick]])) # Unwrap phase if self.last_drift is None: # initial value stored in self.drift_offset self.drift_offset = np.round(self.drift_offset - drift_delta) self.last_drift = self.drift_offset + drift_delta if (self.drift_offset + drift_delta) > (self.last_drift + 0.5): self.drift_offset -= 1 print "(offset -1 sec to %d)" % self.drift_offset elif (self.drift_offset + drift_delta) < (self.last_drift - 0.5): self.drift_offset += 1 print "(offset +1 sec to %d)" % self.drift_offset self.last_drift = drift = self.drift_offset + drift_delta # Make the time be neat; multiple of 3 seconds from midnight # Also don't repeat timestamps; if this t is the same as the last, # drop this one; if it's 6 seconds in between, fill in the gap # This might happen because the chunks are going to beat with the exact 3 seconds t = round_time_to_three_seconds(t) if last_time is not None: if t == last_time: print "------ skipping repeated time %s" % t continue if t == last_time + timedelta(seconds=6): print "------------- filling gap before %s" % t yield {"time": t - timedelta(seconds = 3), "drift": drift, "amplitude": amplitude} last_time = t yield {"time": t, "drift": drift, "amplitude": amplitude} def generate_edge_groups(self, fit_decay=False): """Read samples from source, and generate indices of edge groups""" while True: # Read some samples num_samples = 6 * self.source.fs try: samples = -self.source.get_samples(num_samples) except EOFError: break # Analyse to find edges i_pos_pps, i_neg_pps = self.find_edges(samples[:,self.source.CHANNEL_PPS ]) i_pos_tick, i_neg_tick = self.find_edges(samples[:,self.source.CHANNEL_TICK]) # Fit decay to improve accuracy if required if fit_decay: i_pos_pps = self.fit_decays( samples[:,self.source.CHANNEL_PPS ], i_pos_pps) i_neg_pps = self.fit_decays(-samples[:,self.source.CHANNEL_PPS ], i_neg_pps) i_pos_tick = self.fit_decays( samples[:,self.source.CHANNEL_TICK], i_pos_tick) i_neg_tick = self.fit_decays(-samples[:,self.source.CHANNEL_TICK], i_neg_tick) # Find which is the 'down' tick: the IR signal looks like this: # # Tick: ,u' ,d' ,u' ,d' # # (a): | 0 1 2 3 (long gap) | # (b): | 0 1 (long gap) 2 3 | # # The 3 seconds start either before an 'up' tick (a) or a 'down' tick (b). # (down means pendulum is about to pass through centre, up is return swing) # We distinguish the two cases by comparing the inner gap ([2] - [1]) with # the outer gap ([0]+3secs - [3]). if len(i_pos_tick) < 3: print "Not enough ticks" if PLOT: #plt.clf() fig, ax = plt.subplots(2, sharex=True) ax[0].plot(samples[:, 0]) #[:,self.source.CHANNEL_TICK]) ax[1].plot(samples[:, 1]) ax[0].set_title('tick') ax[1].set_title('pps') for i in i_pos_tick: ax[0].axvline(i, c='g') for i in i_neg_tick: ax[0].axvline(i, c='r') for i in i_pos_pps: ax[1].axvline(i, c='g') for i in i_neg_pps: ax[1].axvline(i, c='r') plt.show() plt.draw() self.source.consume(len(samples)) continue # Start of pulse is up tick_gap_1 = i_pos_tick[1] - i_pos_tick[0] tick_gap_2 = i_pos_tick[2] - i_pos_tick[1] if tick_gap_1 > tick_gap_2: # down tick 0 is the start of the down-swing iref = i_pos_tick[0] print "First tick is down-swing" else: # down tick 1 is the start of the down-swing iref = i_pos_tick[1] print "Second tick is down-swing" # Consume samples belonging to this chunk i_put_back = iref + (3.0 - self.pretrigger)*self.source.fs i_pos_tick = [i for i in i_pos_tick if i >= iref] i_neg_tick = [i for i in i_neg_tick if i >= iref] i_pos_pps = [i for i in i_pos_pps if i >= iref] i_neg_pps = [i for i in i_neg_pps if i >= iref] # XXX This is messy, checking multiple times if len(i_pos_tick) < 3: print "Not enough ticks after down-swing" self.source.consume(len(samples)) continue if PLOT: #plt.clf() fig, ax = plt.subplots(2, sharex=True) ax[0].plot(samples[:, 0]) #[:,self.source.CHANNEL_TICK]) ax[1].plot(samples[:, 1]) ax[0].set_title('tick') ax[1].set_title('pps') for i in i_pos_tick: ax[0].axvline(i, c='g') for i in i_neg_tick: ax[0].axvline(i, c='r') for i in i_pos_pps: ax[1].axvline(i, c='g') for i in i_neg_pps: ax[1].axvline(i, c='r') ax[0].axvline(i_put_back, c='k', lw='2') ax[0].plot(int(iref), samples[int(iref), self.source.CHANNEL_TICK], 'o') ax[0].axvline(iref, ls='--', c='k') plt.show() plt.draw() # Time of reference tick t = self.source.time + timedelta(seconds = iref / self.source.fs) print "Consuming %d samples" % i_put_back #self.put_back_samples(samples[i_put_back:]) self.source.consume(i_put_back) yield t, (i_pos_pps, i_neg_pps), (i_pos_tick, i_neg_tick) def find_edges(self, samples): above = (samples > self.edge_level).astype(int) below = (samples < -self.edge_level).astype(int) i_pos = np.where(np.diff(above) > 0)[0] i_neg = np.where(np.diff(below) > 0)[0] return i_pos, i_neg def fit_decays(self, y, i_edges): """ Fit the exponential decays found in ``y`` after each index in ``ii`` """ Nfit = int(self.decay_fit_duration * self.source.fs) Nlag = int(self.decay_fit_delay * self.source.fs) # Extract decay segments and fit line to log(y) segments = np.vstack([y[i+Nlag:i+Nlag+Nfit] for i in i_edges if (i+Nlag+Nfit) < len(y)]).T K, A_log = np.polyfit(np.arange(Nfit), np.log(segments), 1) A = np.exp(A_log) # Time when fitted line crosses a threshold # when ythresh = A exp(K i) ===> log(ythresh/A)/K = i #i_decay = np.log(0.5 / A) / K i_decay = (np.log(0.5) - A_log) / K return i_edges[:len(i_decay)] + Nlag + i_decay def sanity_check_pps(self, i_edges): # Find mean sample rate from PPS signal fs_mean = np.mean(np.diff(i_edges)) fs_var = np.std(np.diff(i_edges)) print i_edges, fs_mean, fs_var if abs(fs_mean / self.source.fs - 1) > PPS_RATE_ERROR_THRESHOLD: raise DataError("Sample rate is off by too much: %+d ppm" % (1e6*abs(fs_mean/self.source.fs-1))) # PPS should be +/- 1us. Warn if std.deviation * 3, say, is greater than this. # i.e. 9*variance > (1e-6 * fs)^2? but this is rather less than 1 sample... if fs_var > 2: # warn if sample rate seems too variable (empirical) print "** Warning: PPS signal interval variance is high (%d)" % fs_var def calculate_drift(self, ticks, pps, relative_to_pps=False): pps_ref = [i for i in pps if i >= ticks[0]] if not pps_ref: raise DataError("No PPS found after down tick") local_fs = self.source.fs if relative_to_pps: # Assume the gap between PPS edges is 1 second, rather than using # nominal sample rate i_pps_ref = [j for j in range(len(pps)) if pps[j] >= ticks[0]][0] if i_pps_ref > 0: # use previous gap local_fs = pps[i_pps_ref] - pps[i_pps_ref-1] elif (i_pps_ref == 0) and (i_pps_ref+1 < len(pps)): # use next gap local_fs = pps[i_pps_ref+1] - pps[i_pps_ref] # Otherwise fall back on nominal rate drift = (pps_ref[0] - ticks[0]) / local_fs if drift > 1.5: # must be a missing PPS, or something's gone wrong raise DataError("Time between down tick and next PPS too high: %.2f s" % drift) return drift def calculate_amplitude(self, pos, neg): period = (pos[2] - pos[0]) # in samples w = 2 * pi / period # pendulum frequency in rad/sample num = sin(w*pos[0]) - sin(w*neg[0]) + sin(w*pos[1]) - sin(w*neg[1]) den = cos(w*pos[0]) - cos(w*neg[0]) + cos(w*pos[1]) - cos(w*neg[1]) t0 = np.arctan2(num, den) return abs(self.shim_width / (sin(w*pos[0] - t0) - sin(w*neg[0] - t0))) PK[*J clocklogger/input.pyfrom __future__ import division import numpy as np from datetime import datetime, timedelta import pyaudio import logging logger = logging.getLogger(__name__) class PrerecordedDataSource(object): CHANNEL_TICK = 0 CHANNEL_PPS = 1 def __init__(self, filename): self.filename = filename logger.info("Loading pre-recorded data from %s...", filename) data = np.load(filename) self.fs = data['fs'] self.y = data['signal'] self.start_time = datetime.fromtimestamp(data['start_time']) self.i = 0 def get_samples(self, num_samples): """Return some samples""" if self.i >= self.y.shape[0]: raise EOFError samples = self.y[self.i : (self.i + num_samples)] return samples def consume(self, num_samples): """Mark num_samples as having been used""" self.i += num_samples @property def time(self): """Time of first available sample""" return self.start_time + timedelta(seconds = self.i / self.fs) class SoundCardDataSource(object): CHANNEL_TICK = 0 CHANNEL_PPS = 1 def __init__(self, sampling_rate=44100): self.fs = sampling_rate logger.info("Starting PyAudio...") self.pyaudio_manager = pyaudio.PyAudio() dev = self.pyaudio_manager.get_default_input_device_info() if not self.pyaudio_manager.is_format_supported( rate=sampling_rate, input_device=dev['index'], input_channels=2, input_format=pyaudio.paInt16): raise RuntimeError("Unsupported audio format or rate") self.stream = self.pyaudio_manager.open( format=pyaudio.paInt16, channels=2, rate=sampling_rate, input=True) logger.info("PyAudio ready") self.buffer = np.empty((0, 2)) self.buffer_start_time = None def __del__(self): logger.info("Stopping PyAudio stream") self.stream.stop_stream() self.stream.close() self.pyaudio_manager.terminate() def read(self, num_samples): logger.debug("Trying to read %d samples, %d available...", num_samples, self.stream.get_read_available()) raw_data = self.stream.read(num_samples) samples = (np.frombuffer(raw_data, dtype=np.int16) .reshape((-1, 2)) .astype(float) / 2**15) logger.debug("Read %d samples, now %d available", samples.shape[0], self.stream.get_read_available()) return samples def get_samples(self, num_samples): """Return some samples""" num_to_read = num_samples - self.buffer.shape[0] if num_to_read > 0: new_samples = self.read(num_to_read) self.buffer = np.r_[ self.buffer, new_samples ] self.buffer_start_time = \ datetime.utcnow() - timedelta(seconds=self.buffer.shape[0]/self.fs) return self.buffer[:num_samples] def consume(self, num_samples): """Mark num_samples as having been used""" num_to_read = num_samples - self.buffer.shape[0] if num_to_read > 0: self.get_samples(num_to_read) assert self.buffer.shape[0] >= num_samples self.buffer = self.buffer[num_samples:] @property def time(self): """Time of first available sample""" return self.buffer_start_time PKF*J[f@ @ clocklogger/logger.pyimport argparse import time import logging from .input import PrerecordedDataSource, SoundCardDataSource from .analysis import ClockAnalyser, DataError from .output.textfile import TextFileWriter # from .output.influxdb import InfluxDBWriter # from .output.tempodb import TempoDBWriter logger = logging.getLogger(__name__) def get_last_drift(): try: with open('data/last_drift', 'rt') as f: last_drift = float(f.read()) except: last_drift = 0.0 return last_drift def save_last_drift(drift): with open('data/last_drift', 'wt') as f: f.write(str(drift)) def process(analyser, writers): for data in analyser.process(pps_edge='down'): for writer in writers: try: writer.write(data) except Exception as e: logger.error("Writer error [%s]: %s", writer.__class__, e) save_last_drift(data['drift']) def main(): # Set up logging parser = argparse.ArgumentParser(description='clocklogger') parser.add_argument('-L', '--log-level', default='warning') args = parser.parse_args() numeric_level = getattr(logging, args.log_level.upper(), None) if not isinstance(numeric_level, int): raise ValueError('Invalid log level: %s' % args.log_level) logging.basicConfig( level=numeric_level, format="%(asctime)s %(name)s [%(levelname)s] %(message)s") #source = PrerecordedDataSource('../../dataq/record_20130331_0002_100s.npz') source = SoundCardDataSource() analyser = ClockAnalyser(source, initial_drift=get_last_drift()) # Outputs columns = ['time', 'drift', 'amplitude'] # TODO: should do this in a more flexible way writers = [] def add_writer(cls, *args): try: writers.append(cls(*args)) except Exception as err: logger.error("Error creating %s: %s", cls, err) add_writer(TextFileWriter, 'data', 'clock', columns) #add_writer(InfluxDBWriter, 'clock', columns) #add_writer(TempoDBWriter, 'clock', columns) # Read samples, analyze while True: try: process(analyser, writers) except DataError as err: logger.error("Error: %s. Trying to start again in 3 seconds...", err) time.sleep(3) if __name__ == "__main__": main() PK *Jsm--clocklogger/weatherlogger.pyimport time import argparse import logging from datetime import datetime from .source.weather import WeatherStationDataSource from .output.textfile import TextFileWriter from .output.influxdb import InfluxDBWriter from .output.tempodb import TempoDBWriter logger = logging.getLogger(__name__) def round_time_to_interval(t, interval): return t.replace(second=int(t.second//interval)*interval, microsecond=0) def process(source, writers): data = source.get_measurements() data['time'] = round_time_to_interval(datetime.utcnow(), 30) # TODO: interval for writer in writers: try: writer.write(data) except Exception as e: logger.error("Writer error [%s]: %s", writer.__class__, e) def sleep_til_next_time(interval): now = time.time() time.sleep(interval - (now % interval)) def main(): # Set up logging parser = argparse.ArgumentParser(description='weatherlogger') parser.add_argument('-L', '--log-level', default='warning') args = parser.parse_args() numeric_level = getattr(logging, args.log_level.upper(), None) if not isinstance(numeric_level, int): raise ValueError('Invalid log level: %s' % args.log_level) logging.basicConfig( level=numeric_level, format="%(asctime)s %(name)s [%(levelname)s] %(message)s") fields = ['inTemp', 'inHumidity', 'pressure'] source = WeatherStationDataSource(fields) interval = 30.0 # seconds # Outputs columns = ['time'] + fields # TODO: should do this in a more flexible way writers = [] def add_writer(cls, *args): try: writers.append(cls(*args)) except Exception as err: logger.error("Error creating %s: %s", cls, err) add_writer(TextFileWriter, 'data', 'weather', columns) #add_writer(InfluxDBWriter, 'weather', columns) #add_writer(TempoDBWriter, 'weather', columns) # Read data & output while True: process(source, writers) sleep_til_next_time(interval) if __name__ == "__main__": main() PK[*Jclocklogger/output/__init__.pyPK[*JuAAclocklogger/output/influxdb.pyfrom datetime import datetime import logging logger = logging.getLogger(__name__) try: from influxdb import client as influxdb except ImportError: influxdb = None def datetime_to_epoch(d): return int((d - datetime(1970, 1, 1)).total_seconds()) class InfluxDBWriter(object): def __init__(self, base_key, columns): if influxdb is None: raise ImportError("Could not import influxdb package") self.base_key = base_key self.columns = columns self.db = influxdb.InfluxDBClient('localhost', 8086, 'clocklogger', 'pendulum', 'clock') logger.info("Opened connection to InfluxDB") def write(self, data): data = dict(data) data['time'] = datetime_to_epoch(data['time']) points = [ { "name": "clock", "columns": ["%s.%s" % (self.base_key, k) for k in self.columns], "points": [[data[k] for k in self.columns]], } ] self.db.write_points_with_precision(points, 's') PK[*J_+ߘclocklogger/output/tempodb.pyfrom __future__ import absolute_import import os import os.path from tempodb.client import Client from tempodb.protocol import DataPoint import logging logger = logging.getLogger(__name__) class TempoDBWriter(object): DATABASE_ID = "clock" def __init__(self, base_key, columns): try: api_key = os.environ['TEMPODB_API_KEY'] api_sec = os.environ['TEMPODB_API_SECRET'] except KeyError: raise RuntimeError("You must define environment variables " "TEMPODB_API_KEY and TEMPODB_API_SECRET") self.base_key = base_key self.columns = columns self.client = Client(self.DATABASE_ID, api_key, api_sec) def write(self, data): t = data['time'] logger.debug("Data: %s", data) points = [DataPoint.from_data(t, float(data[k]), key='%s.%s' % (self.base_key, k)) for k in self.columns if k != 'time'] resp = self.client.write_multi(points) if resp.status != 200: raise Exception("TempoDB error [%d] %s" % (resp.status, resp.error)) PK[*Jclocklogger/output/textfile.pyimport os import os.path from datetime import datetime import numpy as np import logging logger = logging.getLogger(__name__) def datetime_to_epoch(d): return int((d - datetime(1970, 1, 1)).total_seconds()) class TextFileWriter(object): def __init__(self, path, prefix, columns=None): self.path = path self.columns = columns self.pattern = "%Y/%m/{}%Y-%m-%d.txt".format(prefix) self.file = None def __del__(self): if self.file is not None: self.file.close() def write(self, data): cols = self.columns or sorted(data) formats = { np.float64: '%.6f', float: '%.6f', } time = data['time'] #data['time'] = data['time'].isoformat() data = dict(data) data['time'] = datetime_to_epoch(time) fn = os.path.join(self.path, time.strftime(self.pattern)) if self.file is None or self.file.name != fn: if self.file is not None: self.file.close() if not os.path.exists(os.path.dirname(fn)): os.makedirs(os.path.dirname(fn)) file_already_existed = os.path.exists(fn) self.file = open(fn, 'at') logger.info("Opened %s file %s", "existing" if file_already_existed else "new", fn) self.file.write(" ".join(formats.get(type(data[k]), "%s") % data[k] for k in cols) + "\n") self.file.flush() PK[*Jclocklogger/source/__init__.pyPK[*Jaaclocklogger/source/weather.pyimport logging logger = logging.getLogger(__name__) def _filter_fields(d, fields): """Filter field names in dictionary""" return dict([item for item in d.items() if item[0] in fields]) def _get_driver_class(): try: from weewx.drivers import ws23xx # Patch logging functions to use normal logger not syslog ws23xx.logdbg = lambda msg: logger.debug("WS23xx: %s", msg) ws23xx.loginf = lambda msg: logger.info("WS23xx: %s", msg) ws23xx.logerr = lambda msg: logger.error("WS23xx: %s", msg) ws23xx.logcrt = lambda msg: logger.critical("WS23xx: %s", msg) return ws23xx.WS23xx except ImportError: raise ImportError("WS23xx driver not available") class WeatherStationDataSource(object): def __init__(self, fields=None): """Read weather data from weather station""" # Fill out required config information config = { 'StdArchive': { 'record_generation': 'hardware' } } driver_class = _get_driver_class() self.driver = driver_class(altitude=0.0, config_dict=config) self.fields = fields def get_measurements(self): for packet in self.driver.genLoopPackets(): if self.fields is not None: packet = _filter_fields(packet, self.fields) return packet PK!Hi;h*clocklogger-0.1.dist-info/entry_points.txtN+I/N.,()*OM,H-OOO-RUHOΆPr3󸐤#PK!H|&Ubclocklogger-0.1.dist-info/WHEEL HM K-*ϳR03rOK-J,/RH,Q034 /, (-JLR()*M ILR(4KM̫#DPK!HXX"clocklogger-0.1.dist-info/METADATAM0E@\v%a#Q!FhE[ >6{ϝylct"VxJġtstבe_yDO@ƒda`($KhXZ%[ѐr!r,7eu.Y6^[;0"(\p |@UPLN8|dPK!Hѥ clocklogger-0.1.dist-info/RECORDɒJ}? V3& *( `yh2+B+jիNjj4"79A裻S8 XȒ&!3QEfWR(Z WY&6QƪiJ`D/ZE;a/qáh?^68ׇ'-FGP:'b^~8Hks/I٭&U|/u+CZ颦Mq<oEyva(y6‡gVǭ) vհ)MdI Tn2jbeg3X~Vđ(7_‹fw`Q{ a$8e-Xf.s2V sW!`\L)$ki8ۍ:hBM#&CI1$Œ7U\P"<ʋzx\sbsq|W ˸7\V>M56Upӎ`P{#0,n>Jm?wB[Y-K1qw93mg5WuA0qW3iHG]7|8in}X|x XC,iN}ge<ào6x? FhXwNMψ0x䒩($JQM .%\~b5љV =QЅ;kp:l7)GPK*J44clocklogger/__init__.pyPK[*J$v,,iclocklogger/analysis.pyPK[*J 큂-clocklogger/input.pyPKF*J[f@ @ >;clocklogger/logger.pyPK *Jsm--Dclocklogger/weatherlogger.pyPK[*JMclocklogger/output/__init__.pyPK[*JuAATMclocklogger/output/influxdb.pyPK[*J_+ߘQclocklogger/output/tempodb.pyPK[*JVclocklogger/output/textfile.pyPK[*J\clocklogger/source/__init__.pyPK[*Jaa]clocklogger/source/weather.pyPK!Hi;h*bclocklogger-0.1.dist-info/entry_points.txtPK!H|&Ub'cclocklogger-0.1.dist-info/WHEELPK!HXX"cclocklogger-0.1.dist-info/METADATAPK!Hѥ dclocklogger-0.1.dist-info/RECORDPKbg