Commit d8d748fc authored by Thomas Bruns's avatar Thomas Bruns
Browse files

added seq_multi_fourparsinefit()

parent 3cd64966
......@@ -9,7 +9,7 @@ auxiliary functions related to sine-approximation
import numpy as np
from scipy import linalg as la
from scipy.sparse import coo_matrix
from scipy.sparse import coo_matrix, hstack
from scipy.sparse.linalg import lsqr as sp_lsqr
import matplotlib.pyplot as mp
import warnings
......@@ -766,7 +766,7 @@ def seq_multi_threeparam_Dmatrix(f,t,periods=1):
D = coo_matrix((data,(ri,ci)),shape=(len(t),col)).tocsr()
return np.array(fr), D
def seq_multi_threeparsinefit(f,y,t,periods=1):
def seq_multi_threeparsinefit(f,y,t,periods=1, D_fr=None):
"""
performs a simultanius, sliced three-parameter fit on a multisine signal y
......@@ -791,8 +791,12 @@ def seq_multi_threeparsinefit(f,y,t,periods=1):
Residuals = Difference =(y-y0)
"""
fr, D = seq_multi_threeparam_Dmatrix(f,t,periods) # calculate the design matrix (as sparse matrix)
if D_fr is None:
fr, D = seq_multi_threeparam_Dmatrix(f,t,periods) # calculate the design matrix (as sparse matrix)
else:
D = D_fr[0]
fr = D_fr[1]
abc = sp_lsqr(D,y)
y0 = D.dot(abc[0])
......@@ -805,10 +809,75 @@ def seq_multi_threeparsinefit(f,y,t,periods=1):
f_ab_c = f_ab_c + [0.0,abc[0][-1],0.0]
f_ab_c = np.array(f_ab_c).reshape((len(f_ab_c)//3,3))
print(f_ab_c)
#print(f_ab_c)
return f_ab_c, y0, y-y0 # coefficients, fit signal, residuals
def seq_multi_fourparam_sineFit(f,y,t,periods=1, tol=1.0e-6, n_max=100):
"""
Fit a multi-sinus-signal in slices in one go. With a frequency correction
universal for all frequncies (time-stretch)
Parameters
----------
f : numpy.array of floats
list of frequencies in the signal
t : numpy.array of floats
timestamps of y (typically seconds)
periods : float, optional
the number of periods of each frequncy used for each fit. The default is 1.
tol : float, optional
frequency correction factor (1-tol) when convergence is assumed default 1.0e-6
n_max : uint, optional
maximum number of iterations performed.
Returns
-------
fr, : 1-d numpy.array of floats
frequencies related to slices
D : 2-d numpy.array of floats
Design matrix for the fit
"""
T = t[-1]-t[0]
dw = 1.0
lam=1.0
itr=0
while (np.abs(dw) > tol) and (itr<n_max):
itr += 1
# Designmatrix for three paramsin/cos
fr,D = seq_multi_threeparam_Dmatrix(f, t, periods=periods)
# initial ab
f_ab_c = seq_multi_threeparsinefit(f, y, t, D_fr=(D,fr))[0]
last_col = np.zeros(len(t))
for fi, omi in zip(f,2*np.pi*f):
tau = 1/fi*periods # slice length in seconds
t_sl = np.array_split(t,np.ceil(T/tau)) # array of slices of sample times
sin_cos=np.array([])
for i,ti in enumerate(t_sl): # loop over slices for one frequncy
si = -f_ab_c[i,1]*omi*ti*np.sin(omi * ti)
co = f_ab_c[i,2]*omi*ti*np.cos(omi * ti)
sin_cos = np.hstack((sin_cos,si+co)) # data vector
last_col = last_col+sin_cos
last_col = coo_matrix((last_col,([i for i in range(len(last_col))],[0 for i in range(len(last_col))])),shape=(len(last_col),1)).tocsr()
D = hstack((D, last_col))
abc = sp_lsqr(D,y)[0]
dw = abc[-1]
f = (1-0.9*dw)*f
lam=lam*(1-dw)
if False:
print("abc:")
print(abc)
print("f=")
print(f)
print("cor. = %f" %(lam))
print("took %d iterations" % (itr))
f_ab_c, y0, res = seq_multi_threeparsinefit(f, y, t, D_fr=None)
return f_ab_c, y0, y-y0 # coefficients, fit signal, residuals
def seq_multi_amplitude(f_ab_c):
"""
return the amplitude(s) of a sequentially fit multi-sine signal,
......
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