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

corrected seq_multi_fourparam_... by implementing a stretch factor

parent 1248407b
......@@ -894,7 +894,7 @@ def seq_multi_threeparsinefit(f,y,t,periods=1, D_fr=None, abc0=None, progressive
return f_ab_c, y0, y-y0 # coefficients, fit signal, residuals
def seq_multi_fourparam_Dmatrix(f,t, abc=None, periods=1, progressive=True):
def seq_multi_fourparam_Dmatrix(f,t, delta, abc=None, periods=1, progressive=True):
"""
Fit a multi-sinus-signal in slices in one go.
......@@ -904,6 +904,8 @@ def seq_multi_fourparam_Dmatrix(f,t, abc=None, periods=1, progressive=True):
list of (adapted) frequencies in the signal
t : numpy.array of floats
timestamps of y (typically seconds)
delta: float (scalar)
common stretch factor for the time-scale or shrink for the frequency
abc :numpy.array of floats
parameters from last fit-iteration
periods : float, optional
......@@ -923,7 +925,7 @@ def seq_multi_fourparam_Dmatrix(f,t, abc=None, periods=1, progressive=True):
ci = [] # initialise column indices
ri = [a for a in range(Nt)]*(2*Nf+2) # row indices
data = [0.0]*len(ri) # initialise column indices
data = [0.0]*len(ri) # initialise data
fr = []
col_a =[] # extra part -A*t*sin(om t)
col_b=[] # extra part B*t*cos(om t)
......@@ -936,7 +938,7 @@ def seq_multi_fourparam_Dmatrix(f,t, abc=None, periods=1, progressive=True):
for fi, omi in zip(f,2*np.pi*f):
col_a =[] # extra part -A*t*sin(om t)
col_b=[] # extra part B*t*cos(om t)
om_t = omi*t # omega*t
om_t = omi*t*delta # omega*t
co = np.cos(om_t)
si = np.sin(om_t)
data[f_count*Nt:(f_count+1)*Nt] = co.tolist() # sparse matrix entries
......@@ -950,6 +952,7 @@ def seq_multi_fourparam_Dmatrix(f,t, abc=None, periods=1, progressive=True):
N_samp = int(tau//np.mean(np.diff(t))) # samples/rows per slice
N_slices = len(t)//N_samp # slices for
slic_inds = np.array_split(np.arange(len(t)), N_slices) # r-indices in the slices
if len(slic_inds[-1] < N_samp): # if last slice too short
slic_inds[-2] = np.concatenate([slic_inds[-2],slic_inds[-1]]) # merge last two slices
......@@ -968,7 +971,7 @@ def seq_multi_fourparam_Dmatrix(f,t, abc=None, periods=1, progressive=True):
for i,s in enumerate(slic_inds):
col_a = col_a + [-abc[2*i]]*len(s) # list of A coeffs of cos
col_b = col_b + [abc[2*i+1]]*len(s) # list of B coeffs of sine
last_col = last_col + t*omi*(np.array(col_b)*co + np.array(col_a)*si) # add up 4 param column
last_col = last_col + om_t*(np.array(col_b)*co + np.array(col_a)*si) # add up 4 param column
# else last column is zero as initialised
c_count += 2*len(slic_inds)
......@@ -1027,19 +1030,22 @@ def seq_multi_fourparam_sineFit(f,y,t,periods=1, progressive=True, tol=1.0e-9, n
abc = (f_abc[0:-1,1:3]).flatten() # a,b,a,b,a,b,
abc = np.hstack([abc, f_abc[-2,-2],0.0]) # ababab, c, dw
delta = 1.0 # frequency/time stretch factor
i = 0
while True: # do-while emulation
fr, D = seq_multi_fourparam_Dmatrix(f,t, abc=abc, periods=periods, progressive=progressive) # calculate the design matrix (as sparse matrix)
print("%d. call to solver" % (i+1))
fr, D = seq_multi_fourparam_Dmatrix(f,t,delta, abc=abc, periods=periods, progressive=progressive) # calculate the design matrix (as sparse matrix)
param = sp_lsqr(D,y,x0=abc, atol=1.0e-10, btol=1.0e-10)
#print("abc: %s" % str(param))
f = f*(1+param[0][-1]) # korrekt the frequencies by the factor \Delta\omega
abc = param[0]
print("\delta\Omega: %g" % param[0][-1])
print("f= %s" % str(f))
delta = delta*(1+abc[-1]) # korrekt the timestamps by the factor 1 over \Delta\omega
print("\delta\Omega: %g" % abc[-1])
print("delta= %s" % str(delta))
print("corrected f: %s" % str(delta*f))
if (param[0][-1]/f[0]<tol) :
print("fit arrived at tolerance")
if (np.abs(abc[-1])<tol) :
print("fit arrived at tolerance after %d iterations" % (i+1))
print("last correction: %g" % abc[-1])
break
if (i>n_max):
......@@ -1048,7 +1054,7 @@ def seq_multi_fourparam_sineFit(f,y,t,periods=1, progressive=True, tol=1.0e-9, n
i += 1
y0 = D.dot(param[0])
fr = fr*(1+param[0][-1])
fr = fr*delta
f_ab_c =[]
k=0
......
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