You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 
 
 
 

193 lines
7.7 KiB

  1. #!/usr/bin/python
  2. # Spectral envelope preprocessor.
  3. # Requires two streams as input: the original raw data, and sinefit data.
  4. from nilmdb.utils.printf import *
  5. from nilmdb.utils.time import timestamp_to_human
  6. import nilmtools.filter
  7. import nilmdb.client
  8. from numpy import *
  9. import scipy.fftpack
  10. import scipy.signal
  11. #from matplotlib import pyplot as p
  12. import bisect
  13. from nilmdb.utils.interval import Interval
  14. def main(argv = None):
  15. # Set up argument parser
  16. f = nilmtools.filter.Filter()
  17. parser = f.setup_parser("Spectral Envelope Preprocessor", skip_paths = True)
  18. group = parser.add_argument_group("Prep options")
  19. group.add_argument("-c", "--column", action="store", type=int,
  20. help="Column number (first data column is 1)")
  21. group.add_argument("-n", "--nharm", action="store", type=int, default=4,
  22. help="number of odd harmonics to compute (default 4)")
  23. group.add_argument("-N", "--nshift", action="store", type=int, default=1,
  24. help="number of shifted FFTs per period (default 1)")
  25. exc = group.add_mutually_exclusive_group()
  26. exc.add_argument("-r", "--rotate", action="store", type=float,
  27. help="rotate FFT output by this many degrees (default 0)")
  28. exc.add_argument("-R", "--rotate-rad", action="store", type=float,
  29. help="rotate FFT output by this many radians (default 0)")
  30. group.add_argument("srcpath", action="store",
  31. help="Path of raw input, e.g. /foo/raw")
  32. group.add_argument("sinepath", action="store",
  33. help="Path of sinefit input, e.g. /foo/sinefit")
  34. group.add_argument("destpath", action="store",
  35. help="Path of prep output, e.g. /foo/prep")
  36. # Parse arguments
  37. try:
  38. args = f.parse_args(argv)
  39. except nilmtools.filter.MissingDestination as e:
  40. rec = "float32_%d" % (e.parsed_args.nharm * 2)
  41. print("Source is %s (%s)" % (e.src.path, e.src.layout))
  42. print("Destination %s doesn't exist" % (e.dest.path))
  43. print("You could make it with a command like:")
  44. print(" nilmtool -u %s create %s %s" % (e.dest.url, e.dest.path, rec))
  45. raise SystemExit(1)
  46. if f.dest.layout_count != args.nharm * 2:
  47. print("error: need", args.nharm*2, "columns in destination stream")
  48. raise SystemExit(1)
  49. # Check arguments
  50. if args.column is None or args.column < 1:
  51. parser.error("need a column number >= 1")
  52. if args.nharm < 1 or args.nharm > 32:
  53. parser.error("number of odd harmonics must be 1-32")
  54. if args.nshift < 1:
  55. parser.error("number of shifted FFTs must be >= 1")
  56. if args.rotate is not None:
  57. rotation = args.rotate * 2.0 * pi / 360.0
  58. else:
  59. rotation = args.rotate_rad or 0.0
  60. # Check the sine fit stream
  61. client_sinefit = nilmdb.client.Client(args.url)
  62. sinefit = nilmtools.filter.get_stream_info(client_sinefit, args.sinepath)
  63. if not sinefit:
  64. raise Exception("sinefit data not found")
  65. if sinefit.layout != "float32_3":
  66. raise Exception("sinefit data type is " + sinefit.layout
  67. + "; expected float32_3")
  68. # Check and set metadata in prep stream
  69. f.check_dest_metadata({ "prep_raw_source": f.src.path,
  70. "prep_sinefit_source": sinefit.path,
  71. "prep_column": args.column,
  72. "prep_rotation": repr(rotation),
  73. "prep_nshift": args.nshift })
  74. # Find the intersection of the usual set of intervals we'd filter,
  75. # and the intervals actually present in sinefit data. This is
  76. # what we will process.
  77. filter_int = f.intervals()
  78. sinefit_int = ( Interval(start, end) for (start, end) in
  79. client_sinefit.stream_intervals(
  80. args.sinepath, start = f.start, end = f.end) )
  81. intervals = nilmdb.utils.interval.intersection(filter_int, sinefit_int)
  82. # Run the process (using the helper in the filter module)
  83. f.process_numpy(process, args = (client_sinefit, sinefit.path, args.column,
  84. args.nharm, rotation, args.nshift),
  85. intervals = intervals)
  86. def process(data, interval, args, insert_function, final):
  87. (client, sinefit_path, column, nharm, rotation, nshift) = args
  88. rows = data.shape[0]
  89. data_timestamps = data[:,0]
  90. if rows < 2:
  91. return 0
  92. last_inserted = [nilmdb.utils.time.min_timestamp]
  93. def insert_if_nonoverlapping(data):
  94. """Call insert_function to insert data, but only if this
  95. data doesn't overlap with other data that we inserted."""
  96. if data[0][0] <= last_inserted[0]:
  97. return
  98. last_inserted[0] = data[-1][0]
  99. insert_function(data)
  100. processed = 0
  101. out = zeros((1, nharm * 2 + 1))
  102. # Pull out sinefit data for the entire time range of this block
  103. for sinefit_line in client.stream_extract(sinefit_path,
  104. data[0, 0], data[rows-1, 0]):
  105. def prep_period(t_min, t_max, rot):
  106. """
  107. Compute prep coefficients from time t_min to t_max, which
  108. are the timestamps of the start and end of one period.
  109. Results are rotated by an additional extra_rot before
  110. being inserted into the database. Returns the maximum
  111. index processed, or None if the period couldn't be
  112. processed.
  113. """
  114. # Find the indices of data that correspond to (t_min, t_max)
  115. idx_min = bisect.bisect_left(data_timestamps, t_min)
  116. idx_max = bisect.bisect_left(data_timestamps, t_max)
  117. if idx_min >= idx_max or idx_max >= len(data_timestamps):
  118. return None
  119. # Perform FFT over those indices
  120. N = idx_max - idx_min
  121. d = data[idx_min:idx_max, column]
  122. F = scipy.fftpack.fft(d) * 2.0 / N
  123. # If we wanted more harmonics than the FFT gave us, pad with zeros
  124. if N < (nharm * 2):
  125. F = r_[F, zeros(nharm * 2 - N)]
  126. # Fill output data.
  127. out[0, 0] = round(t_min)
  128. for k in range(nharm):
  129. Fk = F[2 * k + 1] * e**(rot * 1j * (k+1))
  130. out[0, 2 * k + 1] = -imag(Fk) # Pk
  131. out[0, 2 * k + 2] = real(Fk) # Qk
  132. insert_if_nonoverlapping(out)
  133. return idx_max
  134. # Extract sinefit data to get zero crossing timestamps.
  135. # t_min = beginning of period
  136. # t_max = end of period
  137. (t_min, f0, A, C) = [ float(x) for x in sinefit_line.split() ]
  138. t_max = t_min + 1e6 / f0
  139. # Compute prep over shifted windows of the period
  140. # (nshift is typically 1)
  141. for n in range(nshift):
  142. # Compute timestamps and rotations for shifted window
  143. time_shift = n * (t_max - t_min) / nshift
  144. shifted_min = t_min + time_shift
  145. shifted_max = t_max + time_shift
  146. angle_shift = n * 2 * pi / nshift
  147. shifted_rot = rotation - angle_shift
  148. # Run prep computation
  149. idx_max = prep_period(shifted_min, shifted_max, shifted_rot)
  150. if not idx_max:
  151. break
  152. processed = idx_max
  153. # If we processed no data but there's lots in here, pretend we
  154. # processed half of it.
  155. if processed == 0 and rows > 10000:
  156. processed = rows / 2
  157. printf("%s: warning: no periods found; skipping %d rows\n",
  158. timestamp_to_human(data[0][0]), processed)
  159. else:
  160. printf("%s: processed %d of %d rows\n",
  161. timestamp_to_human(data[0][0]), processed, rows)
  162. return processed
  163. if __name__ == "__main__":
  164. main()