Test and development of seq_multi_... functions

SingleSequentialTest.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Sat Dec 12 18:16:02 2020
+
+@author: bruns01
+Sequentieller Multisinusfit mit einer 
+schwachbesetzten Designmatrix 
+"""
+
+import numpy as np
+from scipy.sparse import coo_matrix
+from scipy.sparse.linalg import lsqr
+import matplotlib.pyplot as mp
+
+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
+    
+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, ])
+
+f = f0 * primes/Np0
+T = Np0 /f0
+f_max = np.amax(f)
+
+t = st.sampletimes(2000,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
+
+#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)
+end = datetime.now()
+
+print("Delta=" + str((end-start)))
+