Proper window advancement & insertion for sine fit; comment out plots

src/sinefit.py

@@ -4,7 +4,7 @@ import nilmtools.filter
 import nilmdb.client
 from numpy import *
 from scipy import *
-import pylab as p
+#import pylab as p
 import operator
 
 def main():
@@ -44,15 +44,17 @@ def process(data, interval, args, insert_function, final):
     # Estimate sampling frequency from timestamps
     fs = 1e6 * (rows-1) / (data[-1][0] - data[0][0])
 
-    # Pull out about 4 periods of data at once
-    N = int(4 * fs / f_expected)
+    # Pull out about 3.5 periods of data at once;
+    # we'll expect to match 3 zero crossings in each window
+    N = max(int(3.5 * fs / f_expected), 10)
 
     # If we don't have enough data, don't bother processing it
     if rows < N:
         return 0
 
-    # Process overlapping chunks
+    # Process overlapping windows
     start = 0
+    num_zc = 0
     while start < (rows - N):
         this = data[start:start+N, column]
         t_min = data[start, 0]/1e6
@@ -61,34 +63,50 @@ def process(data, interval, args, insert_function, final):
         # Do 4-parameter sine wave fit
         (A, f0, phi, C) = sfit4(this, fs)
 
-        # Check bounds.  If frequency is too crazy, ignore it
+        # Check bounds.  If frequency is too crazy, ignore this window
         if f0 < (f_expected/2) or f0 > (f_expected*2):
+            print "frequency", f0, "too far from expected value", f_expected
+            start += N
             continue
 
-        p.plot(arange(N), this)
+        #p.plot(arange(N), this)
+        #p.plot(arange(N), A * cos(f0/fs * 2 * pi * arange(N) + phi) + C, 'g')
 
         # Period starts when the argument of cosine is 3*pi/2 degrees,
         # so we're looking for sample number:
-        #     f0/fs * 2 * pi * n + phi = 3 * pi / 2
-        #     f0/fs * 2 * pi * n = (3 * pi / 2 - phi)
         #     n = (3 * pi / 2 - phi) / (f0/fs * 2 * pi)
         zc_n = (3 * pi / 2 - phi) / (f0 / fs * 2 * pi)
         period_n = fs/f0
+
         # Add periods to make N positive
         while zc_n < 0:
             zc_n += period_n
+
+        last_zc = None
         # Mark the zero crossings until we're a half period away
-        # from the end of the window.
+        # from the end of the window
         while zc_n < (N - period_n/2):
-            p.plot(zc_n, C, 'ro')
+            #p.plot(zc_n, C, 'ro')
+            t = t_min + zc_n / fs
+            insert_function([[t * 1e6, f0, A, C]])
+            num_zc += 1
+            last_zc = zc_n
             zc_n += period_n
-        p.plot(zc_n, C, 'bo')
-        p.plot(min(N,zc_n + period_n/2), C, 'go')
 
-        p.show()
-        start += min(N, round(zc_n + period_n/2))
+        # Advance the window one quarter period past the last marked
+        # zero crossing, or advance the window by half its size if we
+        # didn't mark any.
+        if last_zc is not None:
+            advance = min(last_zc + period_n/4, N)
+        else:
+            advance = N/2
+        #p.plot(advance, C, 'go')
+        #p.show()
+
+        start = int(round(start + advance))
 
-    print "processed",start,"rows"
+    # Return the number of rows we've processed
+    print "Marked", num_zc, "zero-crossings in", start, "rows"
     return start
 
 def sfit4(data, fs):