diff --git a/wetb/signal/fix/_rotor_position.py b/wetb/signal/fix/_rotor_position.py
index bb12f6abaf9f9c2eed26ca82fb97bff6d2d17925..efbb6d7da66c61acd73bc08727f57788317e6151 100644
--- a/wetb/signal/fix/_rotor_position.py
+++ b/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
diff --git a/wetb/signal/tests/test_fix.py b/wetb/signal/tests/test_fix.py
index 36c3296c15865d3d44e6e1e6c506a779f7368bc3..0f584284b2651d55b157603f07654b2a504fd794 100644
--- a/wetb/signal/tests/test_fix.py
+++ b/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