no comment

FindPrimesBokeh.py

+#!/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
+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 = []
+    for fs in f_std:
+        try:
+            f1 = (pri[pri/p0<=fs])[-1]/p0  # lower limit
+            r1 = (f1-fs)/fs
+        except:
+            r1 = 0.5
+        f2 = (pri[pri/p0>fs])[0]/p0  # upper limit
+        r2 = (f2-fs)/fs
+
+        r = r1 if (np.abs(r1)<np.abs(r2)) else r2
+        d.append(r)
+        print( f1 if (np.abs(r1)<np.abs(r2)) else f2)
+
+    source.data = dict(x=f_std,y=np.array(d))
+    base.value = """Base = %d """ % p0
+    
+
+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""")
+
+# 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)
+
+curdoc().add_root(column(inputs, plot, width=800))
+curdoc().title = "Sliders"
\ No newline at end of file

SingleSequentialTest.py

@@ -12,143 +12,42 @@ import numpy as np
 from scipy.sparse import coo_matrix
 from scipy.sparse.linalg import lsqr
 import matplotlib.pyplot as mp
+from datetime import datetime
 
 import SineTools as st
 
-def seq_threeparam_Dmatrix(f,t,periods=1):
-    """
-    Fit a multi-sinus-signal in slices in one go.
-
-    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.
-
-    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]
-    col = 0
-    data = np.array([])
-    ci = []
-    ri = []
-    fr = []
-    # Designmatrix for sin/cos
-    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
-        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)
-            cos = np.cos(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
-            
-            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
-
-            Nri = Nri+len(ti)
-
-    # Designmatrix for Bias
-    data = np.hstack((data,np.ones((len(y)))))
-    ci = ci = ci+[col]*len(y)
-    ri = ri + [i for i in range(len(y))]
-    col = col+1
-    # build sparse matrix
-    # print(col)
-    # print(len(data))
-    # print(len(ri))
-    # print(len(ci))
-    # print(len(fr))
-    
-    D = coo_matrix((data,(ri,ci)),shape=(len(y),col)).tocsr()
-    return np.array(fr), D
-
-def seq_threeparam_multisine_fit(f,y,t,periods=1):
-    """
-    performs a simultanius, sliced three-parameter fit on a multisine signal y
-
-    Parameters
-    ----------
-    f : 1-d numpy.array of floats
-        frequencies in the signal
-    y : 1-d numpy.array
-        samples of the signal
-    t : 1-d numpy.array
-        vector of sample times
-    periods : float
-        (fractional) number of periods in one slice The default is 1.
-
-    Returns
-    -------
-    f_ab_c : 2-d numpy.array of floats
-        frequencies, a/b-coefficients and bias related to given frequencies [f1,f1,f1,f2,f2, ...]
-    y0 : 1-d numpy.array of floats
-        samples of the optimal fit
-    resid : 1-d numpy.array of floats
-        Residuals = Difference =(y-y0)
-
-    """
-    fr, D = seq_threeparam_Dmatrix(f,t,periods)  # calculate the design matrix
-    
-    abc = lsqr(D,y)
-    print(abc[0])
-    print(D.shape)
-    print(abc[0].shape)
-    print("###")
-    y0 = D.dot(abc[0])
-    
-    #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]]
-        k=k+2
-    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)
-    return f_ab_c, y0, y-y0 
-
-from datetime import datetime
+from primes import primes
     
 f0 = 0.1
 Np0 = 7.0
 primes = np.array([7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, ])
+primes = primes[1::2]
 
 f = f0 * primes/Np0
-T = Np0 /f0
+T = 3*Np0 /f0
 f_max = np.amax(f)
 
-t = st.sampletimes(2000,T)
+t = st.sampletimes(200,T)
 p = 45.0+180.0*np.array([i for i in range(len(f))])
-m = [0.33]*len(f)
-y = st.multi_waveform_mp(f,m,p,t)+.3
+m = [1.0]*len(f)
+y = st.multi_waveform_mp(f,m,p,t)+.3 + 0.5*(np.random.rand((len(t)))-0.5)
 
-#mp.plot(y,lw=3)
-#mp.show()
+mp.plot(y,lw=3)
+mp.show()
 
 start = datetime.now()
 for i in range(1):
     print(i)
-    abc = st.seq_multi_threeparsinefit(f,y,t,periods=1)
+    abc = st.seq_multi_threeparsinefit(f,y,t,periods=5)
 end = datetime.now()
 
 print("Delta=" + str((end-start)))
-    
+#print(st.seq_multi_amplitude(abc[0]))   
+mag = st.seq_multi_amplitude(abc[0])
+print(np.mean(mag[:,1]))
+print(np.sqrt(np.var(mag[:,1])))
+mp.figure(1)
+mp.clf()
+mp.plot(y,".",ms=1)
+mp.plot(abc[1],lw=1)
+mp.show()