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.
 
 
 
 

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