Browse Source

Port sfit4 function to Python, start writing sinefit filter

tags/nilmtools-1.0
Jim Paris 11 years ago
parent
commit
1db38cc5da
2 changed files with 136 additions and 5 deletions
  1. +5
    -5
      Makefile
  2. +131
    -0
      src/sinefit.py

+ 5
- 5
Makefile View File

@@ -6,11 +6,11 @@ else
endif endif


test: test:
python src/filter.py \
-u "http://localhost:12380/" \
-U "http://127.0.0.1:12380/" \
/lees-compressor/no-leak/raw/1 \
/lees-compressor/no-leak/raw/4
-@make install >/dev/null 2>&1
nilmtool -u "http://localhost/nilmdb" remove /test/foo -s min -e max
python src/sinefit.py \
-u "http://localhost/nilmdb" -c 4 \
/lees-compressor/no-leak/raw/1 /test/foo


test2: test2:
-@nilmtool destroy /lees-compressor/no-leak/raw/4 || true -@nilmtool destroy /lees-compressor/no-leak/raw/4 || true


+ 131
- 0
src/sinefit.py View File

@@ -0,0 +1,131 @@
#!/usr/bin/python

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

def main():
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)')

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

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

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

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

# If we don't have enough data, don't bother processing it
if rows < N:
return 0

# Process overlapping chunks
start = 0
while start < (rows - N):
this = data[start:start+N, column]
(A, f0, phi, C) = sfit4(this, fs)
start += N

print "processed",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

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)
phi = -arctan2(s[1], s[0]) # eqn B.22
C = s[2]

return (A, f0, phi, C)

if __name__ == "__main__":
main()

Loading…
Cancel
Save