Major rewrite:

removed progressive mode in favor of same slice length for all frequencies and an individual bias for each slice.

SineTools.py

@@ -792,57 +792,55 @@ def seq_multi_threeparam_Dmatrix(f,t,periods=1, progressive=True):
           Design matrix for the fit 
 
     """
+    print("## Warning: parameter 'progressive' is ignored in this version of seq_multi_fourparam_Dmatrix")
+
     Nt = len(t)
     Nf = len(f)
 
-    ci = []  # initialise column indices
-    ri = [a for a in range(len(t))]*(2*len(f)+1) # row indices
-    data = [0.0]*len(ri)  # initialise column indices
     fr = []
-    
-    # Designmatrix for cos/sin 
-    c_count = 0 # counter for columns
-    f_count = 0 # counter for frequencies * 2
-    for fi, omi in zip(f,2*np.pi*f):
-        om_t = omi*t        # omega*t
-        data[f_count*Nt:(f_count+1)*Nt] = np.cos(om_t).tolist() # sparse matrix entries
-        data[(f_count+1)*Nt:(f_count+2)*Nt] = np.sin(om_t).tolist() # sparse matrix entries
-        
-        # now split the vectors over the sparse matrix
-        if progressive:
-            tau = 1/fi*(fi//f[0])*periods # approximately same abs. slice length for all fi  
-        else:
-            tau = 1/fi*periods  # slice length in seconds for periods of fi
-        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
-            slic_inds.pop() # remove old last slice
-        
-        fr = fr+[fi]*len(slic_inds)        
-        # for cosine
-        for i,s in enumerate(slic_inds):
-            ci = ci + [c_count+2*i] *len(s)
-        # for sine
-        for i,s in enumerate(slic_inds):
-            ci = ci + [c_count+(2*i+1)] *len(s)
-        
-        c_count += 2*len(slic_inds)                        
-        f_count += 2
-            
-    # Bias part
-    data[2*Nf*Nt:(2*Nf+1)*Nt] = [1.0]*Nt  
-    ci[2*Nf*Nt:(2*Nf+1)*Nt] = [c_count]*Nt
-    c_count += 1
-
+    tau = periods / f[0]
+    N_samp = int(tau // np.mean(np.diff(t)))  # samples(rows) per slice
+    N_slices = len(t) // N_samp  # number of slices
+    slic_inds = np.array_split(np.arange(len(t)), N_slices)  # vectors of r-indices within the slices
+
+    data = np.zeros((Nt*(2*Nf+1)),dtype=np.float)  # initialise matrix values
+    ri = np.zeros((Nt*(2*Nf+1)),dtype=np.uint)  # initialise row indices
+    ci = np.zeros((Nt*(2*Nf+1)),dtype=np.uint)  # initialise colum indices
+
+    r0=0 # row in rectangular matrix
+    c0=0 # column in rectangular matrix
+    ind0 = 0 # index in sparse matrix
+    for slice in slic_inds:
+        Ns = len(slice)
+        data[ind0:ind0+Ns] = 1.0 # bias part for slice i
+        ri[ind0:ind0+Ns] = r0+np.arange(Ns,dtype=np.uint)
+        ci[ind0:ind0+Ns] = c0
+        ind0 += Ns  # proceed Ns entries in sparse matrix
+        fr.extend([0])
+        c0 += 1 # proceed to next column
+        for j,fi in enumerate(f):
+            om_t = 2*np.pi*fi * t[slice] # angular frequency x timestamp
+            # sine matrix entries for freq fi in slice i
+            data[ind0+2*j*Ns:ind0+(2*j+1)*Ns] = np.sin(om_t)
+            ri[ind0+2*j*Ns:ind0+(2*j+1)*Ns] = r0+np.arange(Ns,dtype=np.uint)
+            ci[ind0+2*j*Ns:ind0+(2*j+1)*Ns] = c0+2*j
+            # cosine matrix entries for freq fi in slice i
+            data[ind0+(2*j+1)*Ns:ind0+(2*j+2)*Ns] = np.cos(om_t)
+            ri[ind0+(2*j+1)*Ns:ind0+(2*j+2)*Ns] = r0+np.arange(Ns,dtype=np.uint)
+            ci[ind0+(2*j+1)*Ns:ind0+(2*j+2)*Ns] = c0+2*j+1
+            # build the vector of frequencies
+            fr.extend([fi,fi])
+
+        ind0 += 2*Ns*Nf
+        c0 += 2*Nf # jump 2*Nf columns to next slice-block
+        r0 += Ns   # jump Ns rows to next slice-block
     # build sparse matrix, init as coo, map to csr
-    D = coo_matrix((data,(ri,ci)),shape=(Nt,c_count)).tocsr()
-    # if False:
-    #     print("seq_multi_threeparam_Dmatrix/ D.shape= %s" % str(D.shape))
-    #     mp.spy(D, markersize=5, marker=".",aspect="auto")
-    #     mp.show()
+    print((len(data),len(ri),len(ci)))
+    D = coo_matrix((data,(ri,ci)),shape=(Nt,N_slices*(2*Nf+1))).tocsr()
+    if False:
+        print("seq_multi_threeparam_Dmatrix/ D.shape= %s" % str(D.shape))
+        mp.spy(D, markersize=1, marker=".",aspect="auto")
+        mp.show()
     return np.array(fr), D
 
 
@@ -881,17 +879,14 @@ def seq_multi_threeparsinefit(f,y,t,periods=1, D_fr=None, abc0=None, progressive
     abc = sp_lsqr(D,y,x0=abc0, atol=1.0e-9, btol=1.0e-9)
     y0 = D.dot(abc[0])
     
-    #print(abc)
-    f_ab_c =[]
-    k=0
-    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] # 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(abc[0])
+    N_col = (2*len(f)+1) # length of a row c,a,b,a,b,..,a,b
+    N_slices = int(len(abc[0])/N_col) # number of slices = number of rows
+    f_c_ab = [0]+list(np.array([[fi,fi] for fi in f]).flatten()) + list(abc[0]) # concatenation of arrays
+    f_c_ab = np.array(f_c_ab).reshape((N_slices+1,N_col))
+
     #print(f_ab_c)
-    return f_ab_c, y0, y-y0  # coefficients, fit signal, residuals
+    return f_c_ab, y0, y-y0  # coefficients, fit signal, residuals
 
 
 def seq_multi_fourparam_Dmatrix(f,t, delta, abc=None, periods=1, progressive=True):
@@ -1083,9 +1078,15 @@ def seq_multi_amplitude(f_ab_c):
 
     """
     #print("Test")
-    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))
+    f = f_ab_c[0,1::2]      # vector of frequencies
+    f_vec = []
+    amp = []
+    for i,fi in enumerate(f):
+        a = f_ab_c[1:,1+2*i] # sine coefficients
+        b = f_ab_c[1:,2+2*i] # cosine coefficient
+        f_vec.extend([fi]*len(a))
+        amp.extend(np.abs(1j*a+b))
+    return np.column_stack((np.array(f_vec), np.array(amp)))
 
 def seq_multi_phase(f_ab_c, deg=True):
     """
@@ -1107,9 +1108,15 @@ def seq_multi_phase(f_ab_c, deg=True):
         frequency and associated initial phase. 
 
     """
-    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))
+    f = f_ab_c[0,1::2]      # vector of frequencies
+    f_vec = []
+    amp = []
+    for i,fi in enumerate(f):
+        a = f_ab_c[1:,1+2*i] # cosine coefficients
+        b = f_ab_c[1:,2+2*i] # sine coefficient
+        f_vec.extend([fi]*len(a))
+        amp.extend(np.angle(1j*b+a,deg=deg))
+    return np.column_stack((np.array(f_vec), np.array(amp)))
 
 def seq_multi_bias(f_ab_c):
     """
@@ -1127,7 +1134,8 @@ def seq_multi_bias(f_ab_c):
         Bias of the signal.
 
     """
-    return f_ab_c[-1][1] 
+
+    return f_ab_c[1:,0]
     
 
 def atan_demodulate(samples, fc=0.0, fs=0.0, fc_correct=True, lamb=1.0):