From 65aaac4a798757b8413fabfa4fa45ad85cad8886 Mon Sep 17 00:00:00 2001
From: Thomas Bruns <thomas.bruns@ptb.de>
Date: Fri, 17 Feb 2023 21:38:28 +0100
Subject: [PATCH] added seq_multi_fscale()

---
 SineTools.py | 35 +++++++++++++++++++++++++++++++++++
 1 file changed, 35 insertions(+)

diff --git a/SineTools.py b/SineTools.py
index 412e3d5..8855131 100755
--- a/SineTools.py
+++ b/SineTools.py
@@ -1137,6 +1137,41 @@ def seq_multi_bias(f_ab_c):
 
     return f_ab_c[1:,0]
     
+    
+def seq_multi_fscale(f_ab_c, periods=0):
+     """
+     calculate a frequency correction (scaling) based on linear phase drift
+     in sequential multisine fitting results for phase delay.
+
+     Parameters
+     ----------
+     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.
+
+     Returns
+     -------
+     None.
+
+     """
+     assert periods > 0, "SineTools.seq_multi_fscale: assertion failed, periods > 0 strictly required"
+
+     f_phi = seq_multi_phase(f_ab_c, deg=True)  # frequencies and phases
+     
+     f = f_phi[np.where(f_phi[:,0] == np.amin(f_phi[:,0])),0][0]  # take vector of lowest frequencies
+     phi = f_phi[np.where(f_phi[:,0] == np.amin(f_phi[:,0])),1][0] # take the phases for those
+     
+     t = np.linspace(0.0,len(f),num=len(f),endpoint=False)/f[0]*periods # time incremented per n-periods
+     dphi = np.diff(phi) 
+     dphi[np.where(np.abs(dphi)>=180)] = dphi[np.where(np.abs(dphi)>=180)]-360*np.sign(dphi[np.where(np.abs(dphi)>=180)]) #unwrap
+     phi = np.hstack(([0.0], np.cumsum(dphi)))  # accumulated phase shift per n-periods over time
+     
+     print("t , phi")
+     print(np.vstack((t,phi)).T)
+     slope = np.polyfit(t,phi,1)[0]/f[0]  # slope of the normalized phase over time straight line per period of f_0
+
+     return (1+slope/360.0) # phase was in deg so normalize to period
+
+
 
 def atan_demodulate(samples, fc=0.0, fs=0.0, fc_correct=True, lamb=1.0):
     """
-- 
GitLab