From 23fa2731c2564124cfe9e6ae25b1300eadeb780b Mon Sep 17 00:00:00 2001
From: "Mads M. Pedersen" <mmpe@dtu.dk>
Date: Wed, 3 May 2017 11:10:12 +0200
Subject: [PATCH] bugfix + improved revolution trigger
---
wetb/hawc2/htc_contents.py | 2 +-
wetb/signal/fix/_rotor_position.py | 78 +++++++++++++++++++++++--
wetb/signal/subset_mean.py | 62 +++++++++++++++-----
wetb/signal/tests/test_fix.py | 2 +
wetb/signal/tests/test_subset_mean.py | 17 +++++-
wetb/wind/__init__.py | 1 +
wetb/wind/turbulence/mann_parameters.py | 11 ++--
7 files changed, 145 insertions(+), 28 deletions(-)
diff --git a/wetb/hawc2/htc_contents.py b/wetb/hawc2/htc_contents.py
index 5337613d..93b58f07 100644
--- a/wetb/hawc2/htc_contents.py
+++ b/wetb/hawc2/htc_contents.py
@@ -120,7 +120,7 @@ class HTCContents(object):
self._add_contents(section)
return section
- def add_line(self, name, values, comments):
+ def add_line(self, name, values, comments=""):
line = HTCLine(name, values, comments)
self._add_contents(line)
return line
diff --git a/wetb/signal/fix/_rotor_position.py b/wetb/signal/fix/_rotor_position.py
index 22aac76d..3876204d 100644
--- a/wetb/signal/fix/_rotor_position.py
+++ b/wetb/signal/fix/_rotor_position.py
@@ -5,6 +5,8 @@ Created on 30. mar. 2017
'''
import numpy as np
from wetb.signal.fit._spline_fit import spline_fit
+from wetb.signal.filters._differentiation import differentiation
+
def fix_rotor_position(rotor_position, sample_frq, rotor_speed, fix_dt=None, plt=None):
"""Rotor position fitted with spline
@@ -33,7 +35,7 @@ def fix_rotor_position(rotor_position, sample_frq, rotor_speed, fix_dt=None, plt
from wetb.signal.subset_mean import revolution_trigger
t = np.arange(len(rotor_position))
- indexes = revolution_trigger(rotor_position[:], sample_frq, rotor_speed, max_no_round_diff=4)
+ indexes = revolution_trigger(rotor_position[:], sample_frq, rotor_speed, max_rev_diff=4)
rp = rotor_position[:].copy()
@@ -69,7 +71,9 @@ def fix_rotor_position(rotor_position, sample_frq, rotor_speed, fix_dt=None, plt
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)%360
+ fit = spline_fit(x,y)(t)
+ fit[differentiation(fit, "left")<0]=np.nan
+ return fit%360
def find_fix_dt(rotor_position, sample_frq, rotor_speed, plt=None):
"""Find the optimal fix_dt parameter for fix_rotor_position (function above).
@@ -95,7 +99,8 @@ def find_fix_dt(rotor_position, sample_frq, rotor_speed, plt=None):
"""
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
+ drp = differentiation(fix_rotor_position(rotor_position, sample_frq, rotor_speed, i))
+ rpm_pos = drp%180 / 360 * sample_frq * 60
return np.sum((rpm_pos - rotor_speed)**2)
best = 27
@@ -110,4 +115,69 @@ def find_fix_dt(rotor_position, sample_frq, rotor_speed, plt=None):
plt.show()
return best
-
\ No newline at end of file
+
+def check_rotor_position(rotor_position, sample_frq, rotor_speed, max_rev_diff=1, plt=None):
+ """Rotor position fitted with spline
+
+ 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]
+ 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
+ Fitted rotor position
+ """
+
+ from wetb.signal.subset_mean import revolution_trigger
+
+ t = np.arange(len(rotor_position))/sample_frq
+ indexes = revolution_trigger(rotor_position[:], sample_frq, rotor_speed, max_rev_diff=max_rev_diff)
+
+ rotor_position = rotor_position[:]%360
+
+ rp_fit = fix_rotor_position(rotor_position, sample_frq, rotor_speed, None)
+
+ if plt:
+ #a = (rp.max()-rp.min())/t.max()
+ #plt.plot(t/sample_frq,rp-t*a,label='Continus rotor position (detrended)')
+ f, axarr = plt.subplots(2)
+ print (rp_fit)
+ axarr[0].plot(t, rotor_position)
+ axarr[0].plot(indexes/sample_frq, rotor_position[indexes],'.')
+ axarr[0].plot(t, rp_fit)
+ axarr[1].plot(t, rotor_speed)
+ axarr[1].plot(indexes/sample_frq,60/( differentiation(indexes)/sample_frq))
+ print (t[170:200])
+ drp = differentiation(rp_fit)
+ #drp[(drp<0)]= 0
+ axarr[1].plot(t, drp%180 / 360 * sample_frq * 60, label='fix')
+ plt.legend()
+
+ i1,i2 = indexes[0], indexes[-1]
+ print ("Rev from rotor position", np.sum(np.diff(rotor_position[i1:i2])%180)/360)
+ print ("Rev from rotor speed", np.mean(rotor_speed[i1:i2])*(i2-i1)/60/sample_frq)
+ print ("Rev from fitted rotor position", np.sum(np.diff(rp_fit[i1:i2])%180)/360)
+
+ print ("Mean RPM from rotor speed", np.mean(rotor_speed))
+ print ("Mean RPM from fitted rotor position", np.sum(np.diff(rp_fit[i1:i2])%180)/360 / ((i2-i1)/60/sample_frq))
+
+ spr1 = (np.diff(indexes)/sample_frq)
+
+ #rs1 = 60/( np.diff(indexes)/sample_frq)
+ spr2 = np.array([60/rotor_speed[i1:i2].mean() for i1,i2 in zip(indexes[:-1], indexes[1:])])
+
+ err = spr1-spr2
+ print (err.max())
\ No newline at end of file
diff --git a/wetb/signal/subset_mean.py b/wetb/signal/subset_mean.py
index 0de3ea89..7578c11a 100644
--- a/wetb/signal/subset_mean.py
+++ b/wetb/signal/subset_mean.py
@@ -9,6 +9,7 @@ import unittest
from wetb.utils.geometry import rpm2rads
from _collections import deque
from tables.tests.test_index_backcompat import Indexes2_0TestCase
+from wetb.signal.filters._differentiation import differentiation
def power_mean(power, trigger_indexes, I, rotor_speed, time, air_density=1.225, rotor_speed_mean_samples=1) :
@@ -136,7 +137,7 @@ def cycle_trigger(values, trigger_value=None, step=1, ascending=True, tolerance=
else:
return np.where((values[1:] < trigger_value - tolerance) & (values[:-1] >= trigger_value + tolerance))[0][::step]
-def revolution_trigger(rotor_position, sample_frq, rotor_speed, max_no_round_diff=1):
+def revolution_trigger(rotor_position, sample_frq, rotor_speed, max_rev_diff=1):
"""Returns one index per revolution (minimum rotor position)
Parameters
@@ -155,16 +156,12 @@ def revolution_trigger(rotor_position, sample_frq, rotor_speed, max_no_round_dif
if isinstance(rotor_speed, (float, int)):
rotor_speed = np.ones_like(rotor_position)*rotor_speed
deg_per_sample = rotor_speed*360/60/sample_frq
- sample_per_round = 1/(rotor_speed/60/sample_frq)
- thresshold = deg_per_sample.max()*2
-
- nround_rotor_speed = np.nansum(rotor_speed/60/sample_frq)
-
- mod = [v for v in [5,10,30,60,90] if v>thresshold][0]
+ thresshold = deg_per_sample.max()*3
+ drp = np.diff(rotor_position)%360
+ rp = rotor_position
+
- nround_rotor_position = np.nansum(np.diff(rotor_position)%mod)/360
- #assert abs(nround_rotor_position-nround_rotor_speed)<max_no_round_diff, "No of rounds from rotor_position (%.2f) mismatch with no_rounds from rotor_speed (%.2f)"%(nround_rotor_position, nround_rotor_speed)
- #print (nround_rotor_position, nround_rotor_speed)
+
rp = np.array(rotor_position).copy()
#filter degree increase > thresshold
@@ -176,10 +173,47 @@ def revolution_trigger(rotor_position, sample_frq, rotor_speed, max_no_round_dif
# Best lower is the first lower after upper
best_lower = lower_indexes[np.searchsorted(lower_indexes, upper_indexes)]
upper2lower = best_lower - upper_indexes
- best_lower = best_lower[upper2lower<upper2lower.mean()*2]
- #max_dist_error = max([np.abs((i2-i1)- np.mean(sample_per_round[i1:i2])) for i1,i2 in zip(best_lower[:-1], best_lower[1:])])
- #assert max_dist_error < sample_frq/5, max_dist_error
- return best_lower
+ trigger_indexes = best_lower[upper2lower<upper2lower.mean()*2].tolist()
+
+ if len(trigger_indexes)>1:
+
+
+
+ rpm_rs = np.array([rev.mean() for rev in np.split(rotor_speed, trigger_indexes)[1:-1]])
+ rpm_i = 1/np.diff(trigger_indexes)*60*sample_frq
+ spr_rs = np.array([rev.mean() for rev in np.split(1/rotor_speed*60*sample_frq, trigger_indexes)[1:-1]])
+ spr_i = np.diff(trigger_indexes)
+
+
+ while np.any(spr_rs>spr_i*1.9):
+ i = np.where(spr_rs>spr_i*1.9)[0][0]
+ if np.abs(spr_i[i-1] + spr_i[i] - spr_rs[i])< np.abs(spr_i[i] + spr_i[i+1] - spr_rs[i]):
+ del trigger_indexes[i]
+ else:
+ del trigger_indexes[i+1]
+ spr_i = np.diff(trigger_indexes)
+ spr_rs = np.array([rev.mean() for rev in np.split(1/rotor_speed*60*sample_frq, trigger_indexes)[1:-1]])
+
+
+
+ # if a revolution is too long add trigger
+ if np.any(rpm_rs>rpm_i*2.1):
+ # >1 missing triggers
+ raise NotImplementedError
+ trigger_indexes.extend([np.mean(trigger_indexes[i:i+2]) for i in np.where(rpm_rs>rpm_i*1.9)[0]])
+ trigger_indexes = np.sort(trigger_indexes).astype(np.int)
+
+
+ i1,i2 = trigger_indexes[0], trigger_indexes[-1]
+ nround_rotor_speed = np.nansum(rotor_speed[i1:i2]/60/sample_frq)
+ #mod = [v for v in [5,10,30,60,90] if v>thresshold][0]
+ nround_rotor_position = len(trigger_indexes)-1 #np.nansum(np.diff(rotor_position[i1:i2])%mod)/360
+ if max_rev_diff is not None:
+ diff_pct = abs(nround_rotor_position-nround_rotor_speed)/nround_rotor_position*100
+ assert diff_pct<max_rev_diff, "No of revolution mismatch: rotor_position (%d), rotor_speed (%.1f), diff %.1f%%"%(nround_rotor_position, nround_rotor_speed, diff_pct)
+
+
+ return trigger_indexes
def revolution_trigger_old(values, rpm_dt=None, dmin=5, dmax=10, ):
diff --git a/wetb/signal/tests/test_fix.py b/wetb/signal/tests/test_fix.py
index 86639862..86c9ce47 100644
--- a/wetb/signal/tests/test_fix.py
+++ b/wetb/signal/tests/test_fix.py
@@ -42,6 +42,8 @@ class TestFix(unittest.TestCase):
self.assertEqual(find_fix_dt(ds.azi, sample_frq, ds.Rot_cor), 4)
+
+
if __name__ == "__main__":
#import sys;sys.argv = ['', 'Test.testRotorPositionFix']
diff --git a/wetb/signal/tests/test_subset_mean.py b/wetb/signal/tests/test_subset_mean.py
index 23642ae7..8c4ecd86 100644
--- a/wetb/signal/tests/test_subset_mean.py
+++ b/wetb/signal/tests/test_subset_mean.py
@@ -92,21 +92,32 @@ class TestSubsetMean(unittest.TestCase):
x2 = np.random.randint(0, len(rotor_position),10)
print (list(x2))
else:
- x1 = [447, 854, 595, 804, 847, 488, 412, 199, 675, 766]
+ x1 = [447, 854, 847, 488, 412, 199, 675]
x2 = [92, 647, 821, 422, 33, 159, 369, 99, 157, 464]
rotor_position[x1] += 360
rotor_position[x2] -= 360
rotor_position[90] = 180
+ rotor_position[270:276] = 15
+ rotor_position[598:602] = [360,360,0,0]
+ rotor_position[748:752] = [360,360,0,0]
indexes = revolution_trigger(rotor_position, 20,1/(90/20/60))
- np.testing.assert_array_equal(indexes, [ 91, 180, 270, 360, 450, 540, 630, 720, 810])
if 0:
import matplotlib.pyplot as plt
plt.plot(rotor_position)
plt.plot(indexes, np.zeros_like(indexes),'.')
plt.show()
+ np.testing.assert_array_equal(indexes, [ 91, 180, 270, 360, 450, 540, 630, 720, 810])
+
+ def test_revolution_trigger_max_rev_diff(self):
+ rotor_position = np.arange(0.,360*10,4)
+ rotor_position += np.random.random(len(rotor_position))
+ rotor_position = rotor_position % 360
+ rotor_speed = 1/(90/20/60)
+ revolution_trigger(rotor_position, 20, rotor_speed*1.02, max_rev_diff=3)
+ self.assertRaisesRegex(AssertionError, "No of revolution mismatch", revolution_trigger, rotor_position, 20, rotor_speed*1.02, max_rev_diff=1)
+
-
if __name__ == "__main__":
unittest.main()
diff --git a/wetb/wind/__init__.py b/wetb/wind/__init__.py
index e69de29b..ae17e9a5 100644
--- a/wetb/wind/__init__.py
+++ b/wetb/wind/__init__.py
@@ -0,0 +1 @@
+from .turbulence import *
\ No newline at end of file
diff --git a/wetb/wind/turbulence/mann_parameters.py b/wetb/wind/turbulence/mann_parameters.py
index b6a1d806..b0a94bc7 100644
--- a/wetb/wind/turbulence/mann_parameters.py
+++ b/wetb/wind/turbulence/mann_parameters.py
@@ -254,20 +254,18 @@ def fit_ae(sf, u, L, G, min_bin_count=None, plt=False):
return ae
-def plot_spectra(ae, L, G, k1, uu, vv, ww=None, uw=None, mean_u=1, log10_bin_size=.2, show=True, plt=None):
+def plot_spectra(ae, L, G, k1, uu, vv=None, ww=None, uw=None, mean_u=1, log10_bin_size=.2, plt=None):
# if plt is None:
# import matplotlib.pyplot as plt
_plot_spectra(k1, uu, vv, ww, uw, mean_u, log10_bin_size, plt)
plot_mann_spectra(ae, L, G, "-", mean_u, plt)
- if show:
- plt.show()
-def _plot_spectra(k1, uu, vv, ww=None, uw=None, mean_u=1, log10_bin_size=.2, plt=None):
+def _plot_spectra(k1, uu, vv=None, ww=None, uw=None, mean_u=1, log10_bin_size=.2, plt=None):
if plt is None:
import matplotlib.pyplot as plt
bk1, buu, bvv, bww, buw = logbin_spectra(k1, uu, vv, ww, uw, log10_bin_size)
def plot(xx, label, color, plt):
- plt.semilogx(bk1, bk1 * xx * 10 ** 0 / mean_u ** 2 , '.' + color)
+ plt.semilogx(bk1, bk1 * xx * 10 ** 0 / mean_u ** 2 , '.' + color, label=label)
plot(buu, 'uu', 'r', plt)
if (bvv) is not None:
plot(bvv, 'vv', 'g', plt)
@@ -280,8 +278,9 @@ def plot_mann_spectra(ae, L, G, style='-', u_ref=1, plt=None, spectra=['uu', 'vv
import matplotlib.pyplot as plt
mf = 10 ** (np.linspace(-4, 1, 1000))
muu, mvv, mww, muw = get_mann_model_spectra(ae, L, G, mf)
+ plt.title("ae: %.3f, L: %.2f, G:%.2f"%(ae,L,G))
if 'uu' in spectra: plt.semilogx(mf, mf * muu * 10 ** 0 / u_ref ** 2, 'r' + style)
- if 'vv' in spectra: plt.semilogx(mf, mf * mvv * 10 ** 0 / u_ref ** 2, 'g' + style)
+ if 'vv' in spectra: plt.semilogx(mf, mf * mvv * 10 ** 0 / u_ref ** 2, 'g' + style)
if 'ww' in spectra: plt.semilogx(mf, mf * mww * 10 ** 0 / u_ref ** 2, 'b' + style)
if 'uw' in spectra: plt.semilogx(mf, mf * muw * 10 ** 0 / u_ref ** 2, 'm' + style)
--
GitLab