Compare commits

..

2 Commits

Author SHA1 Message Date
ce0691d6c4 sineefit: Change sfit4 to fit to \sin instead of \cos
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.
2013-04-27 18:12:20 -04:00
4da658e960 sinefit: move initial estimate into the main iteration loop
Just a little less code.  Same results.
2013-04-27 17:50:23 -04:00

View File

@@ -98,12 +98,12 @@ def process(data, interval, args, insert_function, final):
continue 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') #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: # 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 period_n = fs/f0
# Add periods to make N positive # Add periods to make N positive
@@ -149,9 +149,9 @@ def sfit4(data, fs):
Output: Output:
Parameters [A, f0, phi, C] to fit the equation 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: 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 by Jim Paris
(Verified to match sfit4.m) (Verified to match sfit4.m)
@@ -188,12 +188,11 @@ def sfit4(data, fs):
# if something fails with the least squares fit, etc. # if something fails with the least squares fit, etc.
try: try:
# first guess for A0, B0 using 3-parameter fit (step c) # first guess for A0, B0 using 3-parameter fit (step c)
s = zeros(3)
w = 2*pi*f0 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), 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[0] * t * sin(w*t) + s[1] * t * cos(w*t) ] # eqn B.16
s = linalg.lstsq(D, data)[0] # eqn B.18 s = linalg.lstsq(D, data)[0] # eqn B.18
@@ -202,7 +201,7 @@ def sfit4(data, fs):
## Extract results ## Extract results
A = sqrt(s[0]*s[0] + s[1]*s[1]) # eqn B.21 A = sqrt(s[0]*s[0] + s[1]*s[1]) # eqn B.21
f0 = w / (2*pi) 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] C = s[2]
return (A, f0, phi, C) return (A, f0, phi, C)
except Exception as e: except Exception as e: