From 805d8fb24f1b4d0814dd6ce9e33468d2a25aea81 Mon Sep 17 00:00:00 2001 From: Jim Paris Date: Wed, 3 Apr 2013 21:04:03 -0400 Subject: [PATCH] Add prep --- src/prep.py | 126 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 126 insertions(+) create mode 100755 src/prep.py diff --git a/src/prep.py b/src/prep.py new file mode 100755 index 0000000..d1b0a59 --- /dev/null +++ b/src/prep.py @@ -0,0 +1,126 @@ +#!/usr/bin/python + +# Spectral envelope preprocessor. +# Requires two streams as input: the original raw data, and sinefit data. + +import nilmtools.filter +import nilmdb.client +from numpy import * +import scipy.fftpack +import scipy.signal +from matplotlib import pyplot as p +import bisect + +def main(argv = None): + # Set up argument parser + f = nilmtools.filter.Filter() + parser = f.setup_parser("Spectral Envelope Preprocessor", skip_paths = True) + group = parser.add_argument_group("Prep options") + group.add_argument("-c", "--column", action="store", type=int, + help="Column number (first data column is 1)") + group.add_argument("-n", "--nharm", action="store", type=int, default=4, + help="number of odd harmonics to compute") + exc = group.add_mutually_exclusive_group() + exc.add_argument("-r", "--rotate", action="store", type=float, + help="rotate FFT output by this many degrees") + exc.add_argument("-R", "--rotate-rad", action="store", type=float, + help="rotate FFT output by this many radians") + + group.add_argument("srcpath", action="store", + help="Path of raw input, e.g. /foo/raw") + group.add_argument("sinepath", action="store", + help="Path of sinefit input, e.g. /foo/sinefit") + group.add_argument("destpath", action="store", + help="Path of prep output, e.g. /foo/prep") + + # Parse arguments + try: + args = f.parse_args(argv) + except nilmtools.filter.MissingDestination as e: + rec = "float32_%d" % (e.parsed_args.nharm * 2) + print "Source is %s (%s)" % (e.src.path, e.src.layout) + print "Destination %s doesn't exist" % (e.dest.path) + print "You could make it with a command like:" + print " nilmtool -u %s create %s %s" % (e.dest.url, e.dest.path, rec) + raise SystemExit(1) + + # Check arguments + if args.column is None or args.column < 1: + parser.error("need a column number >= 1") + + if args.nharm < 1 or args.nharm > 32: + parser.error("number of odd harmonics must be 1-32") + + if args.rotate is not None: + rotation = args.rotate * 2.0 * pi / 360.0 + else: + rotation = args.rotate_rad or 0.0 + + # Check the sine fit stream + client_sinefit = nilmdb.client.Client(args.url) + sinefit = nilmtools.filter.get_stream_info(client_sinefit, args.sinepath) + if not sinefit: + raise Exception("sinefit data not found") + if sinefit.layout != "float32_3": + raise Exception("sinefit data type is " + sinefit.layout + + "; expected float32_3") + + # Check and set metadata in prep stream + f.check_dest_metadata({ "prep_raw_source": f.src.path, + "prep_sinefit_source": sinefit.path, + "prep_column": args.column }) + + # Run the processing function on all data + f.process_numpy(process, args = (client_sinefit, sinefit.path, args.column, + args.nharm, rotation)) + +def process(data, interval, args, insert_function, final): + (client, sinefit_path, column, nharm, rotation) = args + rows = data.shape[0] + data_timestamps = data[:,0] + + processed = 0 + out = zeros((1, nharm * 2 + 1)) + # Pull out sinefit data for the entire time range of this block + for sinefit_line in client.stream_extract(sinefit_path, + data[0, 0], data[rows-1, 0]): + # Extract sinefit data to get zero crossing timestamps + (t_min, f0, A, C) = [ float(x) for x in sinefit_line.split() ] + t_max = t_min + 1e6 / f0 + + # Find the indices of data that correspond to (t_min, t_max) + idx_min = bisect.bisect_left(data_timestamps, t_min) + idx_max = bisect.bisect_left(data_timestamps, t_max) + if idx_min >= idx_max: + # something's wonky; ignore this period + continue + if idx_max >= len(data_timestamps): + # max is likely past the end of our chunk, so stop + # processing this chunk now. + break + + # Perform FFT over those indices + N = idx_max - idx_min + d = data[idx_min:idx_max, column] + F = scipy.fftpack.fft(d) / N + + # If we wanted more harmonics than we have, pad with zeros + if N < (nharm * 2): + F = r_[F, zeros(nharm * 2 - N)] + + # Fill output data + out[0, 0] = t_min + for k in range(nharm): + Fk = F[2 * k + 1] * e**(rotation * 1j * k) + out[0, 2 * k + 1] = -imag(Fk) # Pk + out[0, 2 * k + 2] = real(Fk) # Qk + + # Insert it and continue + insert_function(out) + processed = idx_max + + print "Processed", processed, "of", rows, "rows" + return processed + +if __name__ == "__main__": + main()