Add median filter

And adjust the period locator accordingly.
Fitting \sin is the same mathematically, it's just conceptually more
straightforward since we're locating zero crossings anyway.

README.txt

@@ -5,10 +5,10 @@ by Jim Paris <jim@jtan.com>
 Prerequisites:
 
   # Runtime and build environments
-  sudo apt-get install python2.7 python2.7-dev python-setuptools
+  sudo apt-get install python2.7 python2.7-dev python-setuptools python-pip
-  sudo apt-get install python-numpy python-scipy python-matplotlib
+  sudo apt-get install python-numpy python-scipy
 
-  nilmdb (1.5.0+)
+  nilmdb (1.6.3+)
 
 Install:
 

setup.py

@@ -79,6 +79,7 @@ setup(name='nilmtools',
               'nilm-copy-wildcard = nilmtools.copy_wildcard:main',
               'nilm-sinefit = nilmtools.sinefit:main',
               'nilm-cleanup = nilmtools.cleanup:main',
+              'nilm-median = nilmtools.median:main',
               ],
           },
       zip_safe = False,

src/filter.py

@@ -236,8 +236,14 @@ class Filter(object):
         metadata = self._client_dest.stream_get_metadata(self.dest.path)
         if not self.force_metadata:
             for key in data:
-                wanted = str(data[key])
+                wanted = data[key]
+                if not isinstance(wanted, basestring):
+                    wanted = str(wanted)
                 val = metadata.get(key, wanted)
+                # Force UTF-8 encoding for comparison and display
+                wanted = wanted.encode('utf-8')
+                val = val.encode('utf-8')
+                key = key.encode('utf-8')
                 if val != wanted and self.dest.rows > 0:
                     m =  "Metadata in destination stream:\n"
                     m += "  %s = %s\n" % (key, val)

src/median.py

@@ -0,0 +1,43 @@
+#!/usr/bin/python
+import nilmtools.filter, scipy.signal
+
+def main(argv = None):
+    f = nilmtools.filter.Filter()
+    parser = f.setup_parser("Median Filter")
+    group = parser.add_argument_group("Median filter options")
+    group.add_argument("-z", "--size", action="store", type=int, default=25,
+                       help = "median filter size (default %(default)s)")
+    group.add_argument("-d", "--difference", action="store_true",
+                       help = "store difference rather than filtered values")
+
+    try:
+        args = f.parse_args(argv)
+    except nilmtools.filter.MissingDestination as e:
+        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, e.src.layout)
+        raise SystemExit(1)
+
+    meta = f.client_src.stream_get_metadata(f.src.path)
+    f.check_dest_metadata({ "median_filter_source": f.src.path,
+                            "median_filter_size": args.size,
+                            "median_filter_difference": repr(args.difference) })
+
+    f.process_numpy(median_filter, args = (args.size, args.difference))
+
+def median_filter(data, interval, args, insert, final):
+    (size, diff) = args
+    (rows, cols) = data.shape
+    for i in range(cols - 1):
+        filtered = scipy.signal.medfilt(data[:, i+1], size)
+        if diff:
+            data[:, i+1] -= filtered
+        else:
+            data[:, i+1] = filtered
+    insert(data)
+    return rows
+
+if __name__ == "__main__":
+    main()

src/prep.py

@@ -80,7 +80,7 @@ def main(argv = None):
     f.check_dest_metadata({ "prep_raw_source": f.src.path,
                             "prep_sinefit_source": sinefit.path,
                             "prep_column": args.column,
-                            "prep_rotation": rotation })
+                            "prep_rotation": repr(rotation) })
 
     # Run the processing function on all data
     f.process_numpy(process, args = (client_sinefit, sinefit.path, args.column,

src/sinefit.py

@@ -98,12 +98,12 @@ def process(data, interval, args, insert_function, final):
             continue
 
         #p.plot(arange(N), this)
-        #p.plot(arange(N), A * cos(f0/fs * 2 * pi * arange(N) + phi) + C, 'g')
+        #p.plot(arange(N), A * sin(f0/fs * 2 * pi * arange(N) + phi) + C, 'g')
 
-        # Period starts when the argument of cosine is 3*pi/2 degrees,
+        # Period starts when the argument of sine is 0 degrees,
         # so we're looking for sample number:
-        #     n = (3 * pi / 2 - phi) / (f0/fs * 2 * pi)
+        #     n = (0 - phi) / (f0/fs * 2 * pi)
-        zc_n = (3 * pi / 2 - phi) / (f0 / fs * 2 * pi)
+        zc_n = (0 - phi) / (f0 / fs * 2 * pi)
         period_n = fs/f0
 
         # Add periods to make N positive
@@ -149,9 +149,9 @@ def sfit4(data, fs):
 
     Output:
       Parameters [A, f0,  phi, C] to fit the equation
-        x[n] = A * cos(f0/fs * 2 * pi * n + phi) + C
+        x[n] = A * sin(f0/fs * 2 * pi * n + phi) + C
       where n is sample number.  Or, as a function of time:
-        x(t) = A * cos(f0 * 2 * pi * t + phi) + C
+        x(t) = A * sin(f0 * 2 * pi * t + phi) + C
 
     by Jim Paris
     (Verified to match sfit4.m)
@@ -188,12 +188,11 @@ def sfit4(data, fs):
     # if something fails with the least squares fit, etc.
     try:
         # first guess for A0, B0 using 3-parameter fit (step c)
+        s = zeros(3)
         w = 2*pi*f0
-        D = c_[cos(w*t), sin(w*t), ones(N)]
-        s = linalg.lstsq(D, data)[0]
 
-        # Now iterate 6 times (step i)
+        # Now iterate 7 times (step b, plus 6 iterations of step i)
-        for idx in range(6):
+        for idx in range(7):
             D = c_[cos(w*t), sin(w*t), ones(N),
                   -s[0] * t * sin(w*t) + s[1] * t * cos(w*t) ] # eqn B.16
             s = linalg.lstsq(D, data)[0] # eqn B.18
@@ -202,7 +201,7 @@ def sfit4(data, fs):
         ## Extract results
         A = sqrt(s[0]*s[0] + s[1]*s[1]) # eqn B.21
         f0 = w / (2*pi)
-        phi = -arctan2(s[1], s[0]) # eqn B.22
+        phi = arctan2(s[0], s[1]) # eqn B.22 (flipped for sin instead of cos)
         C = s[2]
         return (A, f0, phi, C)
     except Exception as e: