Commit 23a37c26 authored by Thomas Bruns's avatar Thomas Bruns
Browse files

merged

parents 2fff5f97 ac01dfcc
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Dec 22 17:15:45 2020
@author: bruns01
"""
import numpy as np
from bokeh.io import curdoc
from bokeh.layouts import column, row
from bokeh.models import ColumnDataSource, Slider, TextInput,Div
from bokeh.plotting import figure
from primes import primes
# Set up callbacks
def update_data(attrname, old, new):
f_std = np.array([1.0, 1.25, 1.6, 2.0, 2.5, 3.2, 4.0, 5.0, 6.3, 8.0 ])
p0 = primes[new]
# pri = primes[(primes<=10*p0) & (primes>p0)]
pri = primes[(primes<=10*p0)]
d = []
p = []
for fs in f_std:
try:
pr1 = (pri[pri/p0<=fs])[-1] # lower limit prime
f1 = pr1/p0 # lower limit frequency
r1 = (f1-fs)/fs
except:
r1 = 0.5
pr2 = (pri[pri/p0>fs])[0] # upper limit prime
f2 = pr2/p0 # upper limit frequency
r2 = (f2-fs)/fs
r = r1 if (np.abs(r1)<np.abs(r2)) else r2
d.append(r)
p.append((pr1 if (np.abs(r1)<np.abs(r2)) else pr2))
#print( f1 if (np.abs(r1)<np.abs(r2)) else f2)
#print( ("%d/%d" % (pr1,p0)) if (np.abs(r1)<np.abs(r2)) else ("%d/%d" % (pr2,p0)))
source.data = dict(x=f_std,y=np.array(d))
base.value = """Base = %d """ % p0
deviation.text = str(np.amax(np.abs(d))) # display maximum deviation
dd = 100*np.amax(np.abs(np.array(d)))
print("%2.2f %%" % dd)
print(p)
f_std = np.array([1.0, 1.25, 1.6, 2.0, 2.5, 3.2, 4.0, 5.0, 6.3, 8.0 ])
p0 = 7
base = TextInput(value="""Base = xx""")
deviation = Div(text="max")
# Set up plot
source = ColumnDataSource(data=dict(x=[], y=[]))
update_data("offset",0,5)
plot = figure(plot_height=400, plot_width=800, title="prime ratio frequencies",
tools="crosshair,pan,reset,save,wheel_zoom",
x_range=[1.0, 10.0], y_range=[-0.1, 0.2])
plot.circle('x', 'y', source=source, radius=.05)
# Set up widgets
offset = Slider(title="offset", value=5, start=0, end=50, step=1)
offset.on_change('value', update_data)
# Set up layouts and add to document
inputs = row(offset, base, deviation)
curdoc().add_root(column(inputs, plot, width=800))
curdoc().title = "Sliders"
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Dec 22 17:15:45 2020
@author: bruns01
"""
import numpy as np
from bokeh.io import curdoc
from bokeh.layouts import column, row
from bokeh.models import ColumnDataSource, Slider, TextInput,Div
from bokeh.plotting import figure
from primes import primes
# Set up callbacks
def update_data(attrname, old, new):
f_std = np.array([1.0, 1.25, 1.6, 2.0, 2.5, 3.2, 4.0, 5.0, 6.3, 8.0 ])
p0 = primes[new]
# pri = primes[(primes<=10*p0) & (primes>p0)]
pri = primes[(primes<=10*p0)]
d = []
p = []
for fs in f_std:
try:
pr1 = (pri[pri/p0<=fs])[-1] # lower limit prime
f1 = pr1/p0 # lower limit frequency
r1 = (f1-fs)/fs
except:
r1 = 0.5
pr2 = (pri[pri/p0>fs])[0] # upper limit prime
f2 = pr2/p0 # upper limit frequency
r2 = (f2-fs)/fs
r = r1 if (np.abs(r1)<np.abs(r2)) else r2
d.append(r)
p.append((pr1 if (np.abs(r1)<np.abs(r2)) else pr2))
#print( f1 if (np.abs(r1)<np.abs(r2)) else f2)
#print( ("%d/%d" % (pr1,p0)) if (np.abs(r1)<np.abs(r2)) else ("%d/%d" % (pr2,p0)))
source.data = dict(x=f_std,y=np.array(d))
base.value = """Base = %d """ % p0
deviation.text = str(np.amax(np.abs(d))) # display maximum deviation
dd = 100*np.amax(np.abs(np.array(d)))
print("%2.2f %%" % dd)
print(p)
f_std = np.array([1.0, 1.25, 1.6, 2.0, 2.5, 3.2, 4.0, 5.0, 6.3, 8.0 ])
p0 = 7
base = TextInput(value="""Base = xx""")
deviation = Div(text="max")
# Set up plot
source = ColumnDataSource(data=dict(x=[], y=[]))
update_data("offset",0,5)
plot = figure(plot_height=400, plot_width=800, title="prime ratio frequencies",
tools="crosshair,pan,reset,save,wheel_zoom",
x_range=[1.0, 10.0], y_range=[-0.1, 0.2])
plot.circle('x', 'y', source=source, radius=.05)
# Set up widgets
offset = Slider(title="offset", value=5, start=0, end=50, step=1)
offset.on_change('value', update_data)
# Set up layouts and add to document
inputs = row(offset, base, deviation)
curdoc().add_root(column(inputs, plot, width=800))
curdoc().title = "Sliders"
......@@ -14,7 +14,6 @@ from scipy.sparse.linalg import lsqr as sp_lsqr
import matplotlib.pyplot as mp
import warnings
def sampletimes(Fs, T): #
"""
generate a t_i vector with \n
......@@ -201,12 +200,11 @@ def seq_threeparsinefit(y, t, f0, periods=1):
RuntimeWarning,
)
Tau = 1.0 / f0 *periods
dt = t[1] - t[0]
dt = np.mean(np.diff(t))
N = int(Tau // dt) ## samples per section
M = int(t.size // N) ## number of sections or periods
abc = np.zeros((M, 3))
for i in range(int(M)):
ti = t[i * N : (i + 1) * N]
yi = y[i * N : (i + 1) * N]
......@@ -214,6 +212,7 @@ def seq_threeparsinefit(y, t, f0, periods=1):
return abc ## matrix of all fit vectors per period
# four parameter sine-fit (with frequency approximation)
def fourparsinefit(y, t, f0, tol=1.0e-7, nmax=1000):
"""
......@@ -303,12 +302,11 @@ def seq_fourparsinefit(y, t, f0, tol=1.0e-7, nmax=1000, periods=1, debug_plot=Fa
RuntimeWarning,
)
Tau = 1.0 / f0 *periods
dt = t[1] - t[0]
dt = np.mean(np.diff(t))
N = int(Tau // dt) ## samples per section
M = int(t.size // N) ## number of sections or periods
abcd = np.zeros((M, 4))
for i in range(M):
ti = t[i * N : (i + 1) * N]
yi = y[i * N : (i + 1) * N]
......@@ -334,7 +332,6 @@ def seq_fourparsinefit(y, t, f0, tol=1.0e-7, nmax=1000, periods=1, debug_plot=Fa
return abcd ## matrix of all fit vectors per period
# fitting a pseudo-random multi-sine signal with 2*Nf+1 parameters
def multi_threeparsinefit(y, t, f0): # f0 vector of frequencies
"""
......@@ -362,7 +359,7 @@ def multi_amplitude(abc): # abc = [a1,b1 , a2,b2, ...,bias]
abc = [a1,b1 , a2,b2, ...,bias] \n
y = sum_i (a_i*sin(2*pi*f0_i*t) + b_i*cos(2*pi*f0_i*t) + c_i
"""
x = abc[:-1][0::2] + 1j * abc[:-1][1::2] # make complex without Bias (last value)
x = abc[:-1][1::2] + 1j * abc[:-1][0::2] # make complex without Bias (last value)
return np.absolute(x)
......@@ -372,7 +369,7 @@ def multi_phase(abc, deg=False): # abc = [bias, a1,b1 , a2,b2, ...]
abc = [a1,b1 , a2,b2, ...,bias] \n
y = sum_i (a_i*sin(2*pi*f0_i*t) + b_i*cos(2*pi*f0_i*t) + c_i
"""
x = abc[:-1][0::2] + 1j * abc[:-1][1::2] # make complex without Bias (last value)
x = abc[:-1][1::2] + 1j * abc[:-1][0::2] # make complex without Bias (last value)
return np.angle(x, deg=deg)
......@@ -431,7 +428,7 @@ def seq_threeparcounterfit(y, t, f0, diff=False):
if diff=True use differentiation to remove carrier (c.f. source)
"""
Tau = 1.0 / f0
dt = t[1] - t[0]
dt = np.mean(np.diff(t))
N = int(np.floor(Tau / dt)) ## samples per section
M = int(np.floor(t.size / N)) ## number of sections or periods
......@@ -706,7 +703,7 @@ def PR_MultiSine(
) # frequency series, sample rate, sample timestamps, waveform
def seq_multi_threeparam_Dmatrix(f,t,periods=1):
def seq_multi_threeparam_Dmatrix(f,t,periods=1, progressive=True):
"""
Fit a multi-sinus-signal in slices in one go.
......@@ -734,29 +731,33 @@ def seq_multi_threeparam_Dmatrix(f,t,periods=1):
ci = []
ri = []
fr = []
# Designmatrix for sin/cos
# Designmatrix for cos/sin
for fi, omi in zip(f,2*np.pi*f):
Nri = 0 # counter for the current row index
tau = 1/fi*periods # slice length in seconds
if progressive:
tau = np.ceil(f[0]/fi*periods)/fi # approximately same abs. slice length for all fi
else:
tau = 1/fi*periods # slice length in seconds for periods of fi
t_sl = np.array_split(t,np.ceil(T/tau)) # array of slices of sample times
fr = fr + [fi]*len(t_sl) # len(t_sl) times frequency fi in Design matrix
for ti in t_sl: # loop over slices fo one frequncy
sin = np.sin(omi * ti)
for ti in t_sl: # loop over slices for one frequency
# cosine part
cos = np.cos(omi * ti)
data = np.hstack((data,sin)) # data vector
data = np.hstack((data,cos)) # data vector
ci = ci+[col]*len(ti) # column index vector
col = col+1
ri = ri+ [Nri+i for i in range(len(ti))] # row index vector
data = np.hstack((data,cos)) # data vector
#sine part
sin = np.sin(omi * ti)
data = np.hstack((data,sin)) # data vector
ci = ci+[col]*len(ti) # column index vector
col = col+1
ri = ri+ [Nri+i for i in range(len(ti))] # row index vector
Nri = Nri+len(ti)
# Designmatrix for Bias
# Bias part
data = np.hstack((data,np.ones((len(t)))))
ci = ci = ci+[col]*len(t)
ri = ri + [i for i in range(len(t))]
......@@ -766,7 +767,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, D_fr=None, abc0=None):
def seq_multi_threeparsinefit(f,y,t,periods=1, D_fr=None, abc0=None, progressive=True):
"""
performs a simultanius, sliced three-parameter fit on a multisine signal y
......@@ -792,10 +793,10 @@ def seq_multi_threeparsinefit(f,y,t,periods=1, D_fr=None, abc0=None):
"""
if D_fr is None:
fr, D = seq_multi_threeparam_Dmatrix(f,t,periods) # calculate the design matrix (as sparse matrix)
fr, D = seq_multi_threeparam_Dmatrix(f,t,periods, progressive=progressive) # calculate the design matrix (as sparse matrix)
else:
D = D_fr[0]
fr = D_fr[1]
D = D_fr[0] # vector of sparse design matrix
fr = D_fr[1] # vector of frequencies
abc = sp_lsqr(D,y,x0=abc0, atol=1.0e-9, btol=1.0e-9)
y0 = D.dot(abc[0])
......@@ -803,10 +804,10 @@ def seq_multi_threeparsinefit(f,y,t,periods=1, D_fr=None, abc0=None):
#print(abc)
f_ab_c =[]
k=0
for fi in fr:
f_ab_c = f_ab_c+ [fi, abc[0][k],abc[0][k+1]]
for fi in fr: # compile a list of frequencies and coefficients
f_ab_c = f_ab_c+ [fi, abc[0][k],abc[0][k+1]] # [f, a, b]...
k=k+2
f_ab_c = f_ab_c + [0.0,abc[0][-1],0.0]
f_ab_c = f_ab_c + [0.0,abc[0][-1],0.0] # add the bias to the list [0,c,0]
f_ab_c = np.array(f_ab_c).reshape((len(f_ab_c)//3,3))
#print(f_ab_c)
......@@ -847,7 +848,7 @@ def seq_multi_fourparam_sineFit(f,y,t,periods=1, tol=1.0e-6, n_max=100):
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)
fr,D = seq_multi_threeparam_Dmatrix(f, t, periods=periods, progressive=False)
# initial ab
f_ab_c = seq_multi_threeparsinefit(f, y, t, D_fr=(D,fr))[0]
last_col = np.zeros(len(t))
......@@ -895,8 +896,8 @@ def seq_multi_amplitude(f_ab_c):
"""
#print("Test")
N = f_ab_c.shape[0]-1
amp = np.array([np.abs(f_ab_c[i,1]+1j*f_ab_c[i][2]) for i in range(f_ab_c.shape[0]-1)]).reshape((N,1))
N = f_ab_c.shape[0]-1 # exclude the bias
amp = np.array([np.abs(1j*f_ab_c[i,1]+f_ab_c[i][2]) for i in range(f_ab_c.shape[0]-1)]).reshape((N,1))
return np.hstack((f_ab_c[:-1,0].reshape((N,1)), amp))
def seq_multi_phase(f_ab_c, deg=True):
......@@ -907,8 +908,9 @@ def seq_multi_phase(f_ab_c, deg=True):
Parameters
----------
f_ab_c : 2-d numpy array of floats (Nx3)
f_ab_c : 2-d numpy array of floats (Nx3)
f,a,b in a row for several rows, as returned by seq_multi_threeparsinefit.
x=a*cos+b*sin+c
deg : Boolean, optional
Flag whether result is in degrees or rad. The default is True (Degrees).
......@@ -918,7 +920,7 @@ def seq_multi_phase(f_ab_c, deg=True):
frequency and associated initial phase.
"""
N = f_ab_c.shape[0]-1
N = f_ab_c.shape[0]-1 # exclude the bias
phase = np.array([np.angle(1j*f_ab_c[i,1]+f_ab_c[i][2],deg=deg) for i in range(f_ab_c.shape[0]-1)]).reshape((N,1))
return np.hstack((f_ab_c[:-1,0].reshape((N,1)), phase))
......
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