nilmtools/src/sinefit.py

#!/usr/bin/python

# Sine wave fitting.  This runs about 5x faster than realtime on raw data.

import nilmtools.filter
import nilmdb.client
from numpy import *
from scipy import *
#import pylab as p
import operator

def main(argv = None):
    f = nilmtools.filter.Filter()
    parser = f.setup_parser("Sine wave fitting")
    group = parser.add_argument_group("Sine fit options")
    group.add_argument('-c', '--column', action='store', type=int,
                       help='Column number (first data column is 1)')
    group.add_argument('-f', '--frequency', action='store', type=float,
                       default=60.0,
                       help='Approximate frequency (default: %(default)s)')
    group.add_argument('-m', '--min-freq', action='store', type=float,
                       help='Minimum valid frequency '
                       '(default: approximate frequency / 2))')
    group.add_argument('-M', '--max-freq', action='store', type=float,
                       help='Maximum valid frequency '
                       '(default: approximate frequency * 2))')

    # Parse arguments
    try:
        args = f.parse_args(argv)
    except nilmtools.filter.MissingDestination as e:
        rec = "float32_3"
        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, rec)
        raise SystemExit(1)

    if args.column is None or args.column < 1:
        parser.error("need a column number >= 1")
    if args.frequency < 0.1:
        parser.error("frequency must be >= 0.1")
    if args.min_freq is None:
        args.min_freq = args.frequency / 2
    if args.max_freq is None:
        args.max_freq = args.frequency * 2
    if (args.min_freq > args.max_freq or
        args.min_freq > args.frequency or
        args.max_freq < args.frequency):
        parser.error("invalid min or max frequency")

    f.check_dest_metadata({ "sinefit_source": f.src.path,
                            "sinefit_column": args.column })
    f.process_numpy(process, args = (args.column, args.frequency,
                                     args.min_freq, args.max_freq))

def process(data, interval, args, insert_function, final):
    (column, f_expected, f_min, f_max) = args
    rows = data.shape[0]

    # Estimate sampling frequency from timestamps
    fs = 1e6 * (rows-1) / (data[-1][0] - data[0][0])

    # 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 windows
    start = 0
    num_zc = 0
    while start < (rows - N):
        this = data[start:start+N, column]
        t_min = data[start, 0]/1e6
        t_max = data[start+N-1, 0]/1e6

        # Do 4-parameter sine wave fit
        (A, f0, phi, C) = sfit4(this, fs)

        # Check bounds.  If frequency is too crazy, ignore this window
        if f0 < f_min or f0 > f_max:
            print "frequency", f0, "outside valid range", f_min, "-", f_max
            start += N
            continue

        #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:
        #     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
        while zc_n < (N - period_n/2):
            #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

        # 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))

    # Return the number of rows we've processed
    print "Marked", num_zc, "zero-crossings in", start, "rows"
    return start

def sfit4(data, fs):
    """(A, f0, phi, C) = sfit4(data, fs)

    Compute 4-parameter (unknown-frequency) least-squares fit to
    sine-wave data, according to IEEE Std 1241-2010 Annex B

    Input:
      data  vector of input samples
      fs    sampling rate (Hz)

    Output:
      Parameters [A, f0,  phi, C] to fit the equation
        x[n] = A * cos(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

    by Jim Paris
    (Verified to match sfit4.m)
    """
    N = len(data)
    t = linspace(0, (N-1) / fs, N)

    ## Estimate frequency using FFT (step b)
    Fc = fft(data)
    F = abs(Fc)
    F[0] = 0   # eliminate DC

    # Find pair of spectral lines with largest amplitude:
    # resulting values are in F(i) and F(i+1)
    i = argmax(F[0:int(N/2)] + F[1:int(N/2+1)])

    # Interpolate FFT to get a better result (from Markus [B37])
    U1 = real(Fc[i])
    U2 = real(Fc[i+1])
    V1 = imag(Fc[i])
    V2 = imag(Fc[i+1])
    n = 2 * pi / N
    ni1 = n * i
    ni2 = n * (i+1)
    K = ((V2-V1)*sin(ni1) + (U2-U1)*cos(ni1)) / (U2-U1)
    Z1 = V1 * (K - cos(ni1)) / sin(ni1) + U1
    Z2 = V2 * (K - cos(ni2)) / sin(ni2) + U2
    i = arccos((Z2*cos(ni2) - Z1*cos(ni1)) / (Z2-Z1)) / n

    # Convert to Hz
    f0 = i * fs / N

    ## Fit it
    # first guess for A0, B0 using 3-parameter fit (step c)
    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)
    for idx in range(6):
        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
        w = w + s[3]	# update frequency estimate

    ## Extract results
    A = sqrt(s[0]*s[0] + s[1]*s[1]) # eqn B.21
    f0 = w / (2*pi)
    try:
        phi = -arctan2(s[1], s[0]) # eqn B.22
    except TypeError:
        # something broke down, just return zeros
        return (0, 0, 0, 0)
    C = s[2]

    return (A, f0, phi, C)

if __name__ == "__main__":
    main()