|
|
@@ -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() |