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.
 
 
 
 

171 lines
6.6 KiB

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