Commit 1b535791 authored by Thomas Bruns's avatar Thomas Bruns
Browse files

correction of four parameter sine fit

parent e4813518
......@@ -9,6 +9,7 @@ auxiliary functions related to sine-approximation
import scipy as sp
from scipy import linalg as la
import matplotlib.pyplot as mp
def sampletimes(Fs, T): #
"""
......@@ -192,7 +193,6 @@ def fourparsinefit(y,t,f0, tol=1.0e-7, nmax=1000):
"""
abcd = threeparsinefit(y,t,f0)
w = 2*sp.pi*f0
err = 1
i = 0
while ((err > tol) and (i<nmax)):
......@@ -205,13 +205,19 @@ def fourparsinefit(y,t,f0, tol=1.0e-7, nmax=1000):
abcd = (la.lstsq(D.transpose(),y))[0]
dw = abcd[3]
w = w+0.9*dw
err = dw/w
i+=1
err = sp.absolute(dw/w)
assert i < nmax, "iteration error"
return sp.hstack((abcd[0:3],w/(2*sp.pi)))
def calc_fourparsine(abcf,t):
"""
return y = abc[0]*sin(2*pi*f0*t) + abc[1]*cos(2*pi*f0*t) + abc[2]
"""
w0 =2*sp.pi*abcf[3]
return abcf[0]*sp.sin(w0*t)+abcf[1]*sp.cos(w0*t)+abcf[2]
"""
from octave ...
......@@ -238,13 +244,15 @@ endfunction
"""
# periodical sinefit at known frequency
def seq_fourparsinefit(y,t,f0, tol=1.0e-7, nmax=1000):
def seq_fourparsinefit(y,t,f0, tol=1.0e-7, nmax=1000, debug_plot=False):
"""
period-wise sine-fit at a known frequency\n
y vector of sample values \n
t vector of sample times\n
f0 known frequency\n
\n
f0 estimate of excitation frequency\n
nmax maximum of iteration to improve f0 \n
debug_plot Flag for plotting the sequential fit for dubugging \n
\n
returns a (n,3)-matrix of coefficient-triplets [[a,b,c], ...]\n
for y = a*sin(2*pi*f0*t) + b*cos(2*pi*f0*t) + c
"""
......@@ -259,6 +267,23 @@ def seq_fourparsinefit(y,t,f0, tol=1.0e-7, nmax=1000):
ti = t[i*N:(i+1)*N]
yi = y[i*N:(i+1)*N]
abcd[i,:] = (fourparsinefit(yi,ti,f0,tol=tol,nmax=nmax))
if debug_plot:
mp.ioff()
fig = mp.figure("seq_fourparsinefit")
fig.clear()
p1 = fig.add_subplot(211)
p2 = fig.add_subplot(212,sharex=p1)
for i in range(M):
p1.plot(t[i*N:(i+1)*N],y[i*N:(i+1)*N],".")
s = calc_fourparsine(abcd[i,:],t[i*N:(i+1)*N]) # fitted data to plot
p1.plot(t[i*N:(i+1)*N],s,"-")
r = y[i*N:(i+1)*N]-s # residuals to plot
p2.plot(t[i*N:(i+1)*N],r,".")
yi = y[i*N:(i+1)*N]
mp.show()
return abcd ## matrix of all fit vectors per period
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment