Improved fit in fix_rotor_position

wetb/signal/fix/_rotor_position.py

@@ -4,8 +4,9 @@ Created on 30. mar. 2017
 @author: mmpe
 '''
 import numpy as np
-def fix_rotor_position(rotor_position, sample_frq, rotor_speed, polynomial_sample_length=500):
-    """Rotor position fitted with multiple polynomials
+from wetb.signal.fit._spline_fit import spline_fit
+def fix_rotor_position(rotor_position, sample_frq, rotor_speed, fix_dt=None, plt=None):
+    """Rotor position fitted with spline
     
     Parameters
     ----------
@@ -15,9 +16,14 @@ def fix_rotor_position(rotor_position, sample_frq, rotor_speed, polynomial_sampl
         Sample frequency [Hz]
     rotor_speed : array_like
         Rotor speed [RPM]
-    polynomial_sample_length : int, optional
-        Sample lengths of polynomial used for fit
-        
+    fix_dt : int, float or None, optional
+        Time distance [s] between spline fix points\n
+        If None (default) a range of seconds is tested and the result that minimize the RMS 
+        between differentiated rotor position fit and rotor speed is used.\n 
+        Note that a significant speed up is achievable by specifying the parameter
+    plt : PyPlot or None
+        If PyPlot a visual interpretation is plotted
+            
     Returns
     -------
     y : nd_array
@@ -27,48 +33,79 @@ def fix_rotor_position(rotor_position, sample_frq, rotor_speed, polynomial_sampl
     from wetb.signal.subset_mean import revolution_trigger
     
     t = np.arange(len(rotor_position))
-    indexes = revolution_trigger(rotor_position[:].copy(), sample_frq, rotor_speed, max_no_round_diff=4)
+    indexes = revolution_trigger(rotor_position[:], sample_frq, rotor_speed, max_no_round_diff=4)
     
     rp = rotor_position[:].copy()
     
     for i in indexes:
         rp[i:] += 360
         
-    N = polynomial_sample_length
-    N2 = N/2
-    
-    rp_fita = np.empty_like(rp)+np.nan
-    rp_fitb = np.empty_like(rp)+np.nan
+    if fix_dt is None:
+        fix_dt = find_fix_dt(rotor_position, sample_frq, rotor_speed)
+        
+    N = int(np.round(fix_dt*sample_frq))
+    N2 = N//2
+    if plt:
+        a = (rp.max()-rp.min())/t.max()
+        plt.plot(t/sample_frq,rp-t*a,label='Continus rotor position (detrended)')
     
-    for i in range(0,len(rp)+N, N):
+    points = []
+    for j, i in enumerate(range(0,len(rp), N)):
         #indexes for subsets for overlapping a and b polynomials
-        ia1 = max(0,i-N2)
-        ia2 = min(len(rp),i+N2)
-        ib1 = i
-        ib2 = i+N
+        i1 = max(0,i-N2)
+        i2 = min(len(rp),i+N2)
         
         #fit a polynomial
-        if ia1<len(rp):
-            z = np.polyfit(t[ia1:ia2]-t[ia1], rp[ia1:ia2], 3)
-            rp_fita[ia1:ia2] = np.poly1d(z)(t[ia1:ia2]-t[ia1])
-            
-        if ib1<len(rp):
-            #fit b polynomial
-            z = np.polyfit(t[ib1:ib2]-t[ib1], rp[ib1:ib2], 3)
-            rp_fitb[ib1:ib2] = np.poly1d(z)(t[ib1:ib2]-t[ib1])
+        if i1<len(rp):
+            poly_coef = np.polyfit(t[i1:i2]-t[i1], rp[i1:i2], 1)
+            points.append((t[i],np.poly1d(poly_coef)(t[i]-t[i1])))
+            if plt:
+                plt.plot(t[i1:i2]/sample_frq, np.poly1d(poly_coef)(t[i1:i2]-t[i1])- t[i1:i2]*a, 'mc'[j%2], label=('', "Line fit for points (detrended)")[j<2])
+                
+    x,y = np.array(points).T
     
-    weighta = (np.cos(np.arange(len(rp))/N*2*np.pi))/2+.5
-    weightb = (-np.cos(np.arange(len(rp))/N*2*np.pi))/2+.5
+    if plt:
+        plt.plot(x/sample_frq,y-x*a,'.', label='Fit points (detrended)')
+        plt.plot(t/sample_frq, spline_fit(x,y)(t)-t*a, label='Spline (detrended)')
+        plt.legend()
+        plt.show()
+    return spline_fit(x,y)(t)
+
+def find_fix_dt(rotor_position, sample_frq, rotor_speed, plt=None):
+    """Find the optimal fix_dt parameter for fix_rotor_position (function above).
+    Optimal is defined as the value that minimizes the sum of squared differences 
+    between differentiated rotor position and rotor speed 
     
-    return (rp_fita*weighta + rp_fitb*weightb)%360
+    Parameters
+    ----------
+    rotor_position : array_like
+        Rotor position [deg] (0-360)
+    sample_frq : int or float
+        Sample frequency [Hz]
+    rotor_speed : array_like
+        Rotor speed [RPM]
+    plt : pyplot or None
+        If pyplot, a visual interpretation is plotted
+        
+    Returns
+    -------
+    y : int
+        Optimal value for the fix_dt parameter for fix_rotor_position
     
-
-def find_polynomial_sample_length(rotor_position, sample_frq, rotor_speed):
+    """
     from wetb.signal.filters import differentiation
     def err(i):
         rpm_pos = differentiation(fix_rotor_position(rotor_position, sample_frq, rotor_speed, i))%180 / 360 * sample_frq * 60
         return np.sum((rpm_pos - rotor_speed)**2)
-    x_lst = np.arange(100,2000,100)
-    res = [err(x) for x in x_lst]
-    return x_lst[np.argmin(res)]
+    
+    best = 0
+    for r in [np.arange(7,46,12), np.arange(-6,7,3),np.arange(-2,3)]:
+        x_lst = r + best
+        res = [err(x) for x in x_lst]
+        if plt is not None:
+            plt.plot(x_lst, res,'.-')
+        best = x_lst[np.argmin(res)]
+    if plt is not None:
+        plt.show()
+    return best 
     
\ No newline at end of file

wetb/signal/tests/test_fix.py

@@ -6,8 +6,7 @@ Created on 30. mar. 2017
 import os
 import unittest
 from wetb import gtsdf
-from wetb.signal.fix._rotor_position import fix_rotor_position,\
-    find_polynomial_sample_length
+from wetb.signal.fix._rotor_position import fix_rotor_position, find_fix_dt
 from wetb.signal.filters._differentiation import differentiation
 
 
@@ -19,22 +18,26 @@ class TestFix(unittest.TestCase):
     def testRotorPositionFix(self):
         ds = gtsdf.Dataset(tfp+'azi.hdf5')
         sample_frq = 25
-        print (find_polynomial_sample_length(ds.azi, sample_frq, ds.Rot_cor))
-        rp_fit = fix_rotor_position(ds.azi, sample_frq, ds.Rot_cor, 500)
-        
+         
+        import matplotlib.pyplot as plt
+        print (find_fix_dt(ds.azi, sample_frq, ds.Rot_cor, plt))
+        rp_fit = fix_rotor_position(ds.azi, sample_frq, ds.Rot_cor)
+                  
         rpm_pos = differentiation(rp_fit)%180 / 360 * sample_frq * 60
-        
-        self.assertLess(np.sum((rpm_pos - ds.Rot_cor)**2),50)
+        err_sum = np.sum((rpm_pos - ds.Rot_cor)**2)
+         
+        self.assertLess(err_sum,40)
         if 0:
             import matplotlib.pyplot as plt
-            plt.plot( differentiation(ds.azi)%180 / 360 * sample_frq * 60, label='fit')
-            plt.plot(ds.Rot_cor)
-            plt.plot( differentiation(rp_fit)%180 / 360 * sample_frq * 60, label='fit')
+            t = ds.Time-ds.Time[0]
+            plt.plot(t, differentiation(ds.azi)%180 / 360 * sample_frq * 60, label='fit')
+            plt.plot(t, ds.Rot_cor)
+            plt.plot(t, differentiation(rp_fit)%180 / 360 * sample_frq * 60, label='fit')
             plt.ylim(10,16)
-            
+              
             plt.show()
 
-
+ 
 if __name__ == "__main__":
     #import sys;sys.argv = ['', 'Test.testRotorPositionFix']
     unittest.main()
\ No newline at end of file