diff --git a/wetb/dlc/high_level.py b/wetb/dlc/high_level.py
index 6c29c4589245b066667d89bf4b55cf05593eaa87..418de910337dafd4fb557fdb99e75e126d8b1f27 100644
--- a/wetb/dlc/high_level.py
+++ b/wetb/dlc/high_level.py
@@ -44,6 +44,36 @@ def Weibull2(u, k, wsp_lst):
edges = np.r_[wsp_lst[0] - (wsp_lst[1] - wsp_lst[0]) / 2, (wsp_lst[1:] + wsp_lst[:-1]) / 2, wsp_lst[-1] + (wsp_lst[-1] - wsp_lst[-2]) / 2]
return [-cdf(e1) + cdf(e2) for wsp, e1, e2 in zip(wsp_lst, edges[:-1], edges[1:])]
+def Weibull_IEC(Vref, Vhub_lst):
+ """Weibull distribution according to IEC 61400-1:2005, page 24
+
+ Parameters
+ ----------
+ Vref : int or float
+ Vref of wind turbine class
+ Vhub_lst : array_like
+ Wind speed at hub height. Must be equally spaced.
+
+ Returns
+ -------
+ nd_array : list of probabilities
+
+ Examples
+ --------
+ >>> Weibull_IEC(50, [4,6,8])
+ [ 0.11002961 0.14116891 0.15124155]
+ """
+ Vhub_lst = np.array(Vhub_lst)
+ #Average wind speed
+ Vave=.2*Vref
+ #Rayleigh distribution
+ Pr = lambda x : 1 - np.exp(-np.pi*(x/(2*Vave))**2)
+ #Wsp bin edges: [4,6,8] -> [3,5,7,9]
+ wsp_bin_edges = np.r_[Vhub_lst[0] - (Vhub_lst[1] - Vhub_lst[0]) / 2, (Vhub_lst[1:] + Vhub_lst[:-1]) / 2, Vhub_lst[-1] + (Vhub_lst[-1] - Vhub_lst[-2]) / 2]
+ #probabilities of 3-5, 5-7, 7-9
+ return np.array([-Pr(e1) + Pr(e2) for e1, e2 in zip(wsp_bin_edges[:-1], wsp_bin_edges[1:])])
+
+
class DLCHighLevel(object):
@@ -159,7 +189,7 @@ class DLCHighLevel(object):
dist = self.dlc_df[dist_key][row]
if str(dist).lower() == "weibull" or str(dist).lower() == "rayleigh":
- dist = Weibull2(self.vref * .2, self.shape_k, values)
+ dist = Weibull_IEC(self.vref, values)
else:
def fmt(v):
if "#" in str(v):
diff --git a/wetb/dlc/tests/test_high_level.py b/wetb/dlc/tests/test_high_level.py
index 82a658835e669e305807d92ab36f62d26de7f78f..473e81bd6d36eb1adc54304b2d75ede8b0413fdc 100644
--- a/wetb/dlc/tests/test_high_level.py
+++ b/wetb/dlc/tests/test_high_level.py
@@ -10,7 +10,7 @@ from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
import unittest
-from wetb.dlc.high_level import DLCHighLevel, Weibull
+from wetb.dlc.high_level import DLCHighLevel, Weibull, Weibull_IEC
import os
import numpy as np
@@ -109,7 +109,10 @@ class TestDLCHighLevel(unittest.TestCase):
p_tot = np.array([value for key, value in weibull.items()]).sum()
self.assertTrue(np.allclose(p_tot, 1.0))
-
+ def test_weibull_IEC(self):
+ Vref = 50
+ np.testing.assert_array_almost_equal(Weibull_IEC(Vref, [4,6,8]), [ 0.11002961, 0.14116891, 0.15124155])
+
if __name__ == "__main__":
#import sys;sys.argv = ['', 'Test.testName']
diff --git a/wetb/signal/__init__.py b/wetb/signal/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..42d8da654fe6b7fc28c9998298c8d8ef5c397ac5
--- /dev/null
+++ b/wetb/signal/__init__.py
@@ -0,0 +1,12 @@
+
+d = None
+d = dir()
+
+from .interpolation import interpolate
+
+
+__all__ = [m for m in set(dir()) - set(d)]
+
+
+
+
diff --git a/wetb/signal/filters/__init__.py b/wetb/signal/filters/__init__.py
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..ebe5b72cdfc049ea8a1f5a485d5949ccd8782bcd 100644
--- a/wetb/signal/filters/__init__.py
+++ b/wetb/signal/filters/__init__.py
@@ -0,0 +1,2 @@
+from ._despike import *
+from ._differentiation import *
\ No newline at end of file
diff --git a/wetb/signal/filters/despike.py b/wetb/signal/filters/_despike.py
similarity index 85%
rename from wetb/signal/filters/despike.py
rename to wetb/signal/filters/_despike.py
index f7256c68b1bd499faf78987afa83b644ff8f4dee..7fda9e33462d907b229b3339c87a98d0197d0359 100644
--- a/wetb/signal/filters/despike.py
+++ b/wetb/signal/filters/_despike.py
@@ -4,8 +4,7 @@ Created on 13/07/2016
@author: MMPE
'''
import numpy as np
-from wetb.signal.filters.first_order import low_pass
-from wetb.signal.filters import replacer
+from wetb.signal.filters import replacer, frq_filters
replace_by_nan = replacer.replace_by_nan
@@ -42,15 +41,15 @@ def univeral_thresshold_finder(data, variation='mad', plt=None):
thresshold = np.sqrt(2 * np.log(data.shape[0])) * variation # Universal thresshold (expected maximum of n random variables)
return thresshold_finder(data, thresshold, plt)
-def despike(data, dt, spike_finder=univeral_thresshold_finder, spike_replacer=replace_by_nan, plt=None):
+def despike(data, spike_length, spike_finder=univeral_thresshold_finder, spike_replacer=replace_by_nan, plt=None):
"""Despike data
Parameters
---------
data : array_like
data
- dt : int or float
- time step
+ spike_length : int
+ Typical spike duration [samples]
spike_finder : function
Function returning indexes of the spikes
spike_replacer : function
@@ -63,8 +62,10 @@ def despike(data, dt, spike_finder=univeral_thresshold_finder, spike_replacer=re
if plt:
plt.plot(data, label='Input')
data = np.array(data).copy()
- lp_data = low_pass(data, dt, 1)
+ lp_data = low_pass(data, spike_length, 1)
hp_data = data - lp_data
+ hp_data = frq_filters.high_pass(data, spike_length*4, 1, order=2)
+ #spike_length, sample_frq/2, 1, order=1)
spike_mask = spike_finder(hp_data, plt=plt)
despike_data = spike_replacer(data, spike_mask)
if plt:
diff --git a/wetb/signal/filters/_differentiation.py b/wetb/signal/filters/_differentiation.py
new file mode 100644
index 0000000000000000000000000000000000000000..16caa027209d4b5ec1ba3ae59ffcc8be93f14910
--- /dev/null
+++ b/wetb/signal/filters/_differentiation.py
@@ -0,0 +1,37 @@
+'''
+Created on 29. mar. 2017
+
+@author: mmpe
+'''
+from wetb.signal.filters.frq_filters import low_pass
+import numpy as np
+def differentiation(x,sample_frq=None, cutoff_frq=None):
+ """Differentiate the signal
+
+ Parameters
+ ----------
+ x : array_like
+ The input signal
+ sample_frq : int, float or None, optional
+ sample frequency of signal (only required if low pass filer is applied)
+ cutoff_frq : int, float or None, optional
+ Low pass filter cut off (frequencies higher than this frequency will be suppressed)
+ Returns
+ -------
+ y : ndarray
+ differentiated signal
+
+
+ Examples
+ --------
+ >>> differentiation([1,2,1,0,1,1])
+ """
+
+
+ if cutoff_frq is not None:
+ assert sample_frq is not None, "Argument sample_frq must be set to apply low pass filter"
+ x = low_pass(x, sample_frq, cutoff_frq)
+ else:
+ x = np.array(x)
+ dy = np.r_[x[1]-x[0], (x[2:]-x[:-2])/2, x[-1]-x[-2]]
+ return dy
\ No newline at end of file
diff --git a/wetb/signal/filters/frq_filters.py b/wetb/signal/filters/frq_filters.py
new file mode 100644
index 0000000000000000000000000000000000000000..58dec73c5078a9a3bb58709b0cedfd0be500f4a5
--- /dev/null
+++ b/wetb/signal/filters/frq_filters.py
@@ -0,0 +1,144 @@
+'''
+Created on 27. mar. 2017
+
+@author: mmpe
+'''
+import numpy as np
+from scipy import signal
+
+def sine_generator(sample_frq, sine_frq, duration):
+ """Create a sine signal for filter test
+
+ Parameters
+ ----------
+ sample_frq : int, float
+ Sample frequency of returned signal [Hz]
+ sine_frq : int, float
+ Frequency of sine signal [Hz]
+ duration : int, float
+ Duration of returned signal [s]
+
+ Returns
+ -------
+ t,sig : ndarray, ndarray
+ time (t) and sine signal (sig)
+
+ Examples
+ --------
+ >>> sine_generator(10,1,2)
+ """
+ T = duration
+ nsamples = sample_frq * T
+ w = 2. * np.pi * sine_frq
+ t = np.linspace(0, T, nsamples, endpoint=False)
+ sig = np.sin(w * t)
+ return t, sig
+
+
+def low_pass(x, sample_frq, cutoff_frq, order=5):
+ """Low pass filter (butterworth)
+
+ Parameters
+ ----------
+ x : array_like
+ Input signal
+ sample_frq : int, float
+ Sample frequency [Hz]
+ cutoff_frq : int, float
+ Cut off frequency [Hz]
+ order : int
+ Order of the filter (1th order: 20db per decade, 2th order:)
+
+ Returns
+ -------
+ y : ndarray
+ Low pass filtered signal
+
+
+ Examples
+ --------
+ >>>
+ """
+ nyquist_frq = 0.5 * sample_frq
+ normal_cutoff = cutoff_frq / nyquist_frq
+ b,a = signal.butter(order, normal_cutoff, btype='low', analog=False)
+ return signal.filtfilt(b, a, x)
+
+def high_pass(x, sample_frq, cutoff_frq, order=5):
+ """Low pass filter (butterworth)
+
+ Parameters
+ ----------
+ x : array_like
+ Input signal
+ sample_frq : int, float
+ Sample frequency [Hz]
+ cutoff_frq : int, float
+ Cut off frequency [Hz]
+ order : int
+ Order of the filter (1th order: 20db per decade, 2th order:)
+
+ Returns
+ -------
+ y : ndarray
+ Low pass filtered signal
+
+
+ Examples
+ --------
+ >>>
+ """
+ nyquist_frq = 0.5 * sample_frq
+ normal_cutoff = cutoff_frq / nyquist_frq
+ b,a = signal.butter(order, normal_cutoff, btype='high', analog=False)
+ return signal.filtfilt(b, a, x)
+
+def frequency_response(sample_frq, cutoff_frq, type, order, plt=None):
+ """Frequency response of low/high pass filter (butterworth
+
+ Parameters
+ ----------
+ sample_frq : int, float
+ Sample frequency [Hz]
+ cutoff_frq : int, float
+ Cut off frequency [Hz]
+ type : {'low','high'}
+ Low or high pass filter
+ order : int
+ Order of the filter (1th order: 20db per decade, 2th order: 40db per decade)
+ plt : pyplot, optional
+ If specified, the frequency response is plotted
+
+ Returns
+ -------
+ w,h : ndarray, ndarray
+ Frequency (w) in Hz and filter response in db
+
+
+ Examples
+ --------
+ >>>
+ """
+ nyquist_frq = 0.5 * sample_frq
+ normal_cutoff = cutoff_frq / nyquist_frq
+ assert 0<normal_cutoff<1, "cutoff frequency must be <= nyquist frequency"
+ b,a = signal.butter(order, cutoff_frq, btype=type, analog=True)
+ w, h = signal.freqs(b, a)
+ h_db = 20 * np.log10(abs(h))
+ if plt:
+ plt.plot(w, h_db, label='%d order filter response'%order)
+
+ plt.legend(loc=0)
+
+ title = 'Butterworth filter frequency response'
+ if plt.axes().get_title()!=title:
+ plt.title(title)
+ plt.xlabel('Frequency [Hz]')
+ plt.ylabel('Amplitude [dB]')
+ plt.margins(.1, .1)
+ plt.xscale('log')
+
+ plt.grid(which='both', axis='both')
+ plt.axvline(cutoff_frq, color='green') # cutoff frequency
+
+ return w,h_db
\ No newline at end of file
diff --git a/wetb/signal/fit/__init__.py b/wetb/signal/fit/__init__.py
index 2578ab2dde2c387e6801224980c9c13848c1b548..4e0ddd1f8dc62eccddc5845fe9e21e9f38cdc892 100644
--- a/wetb/signal/fit/__init__.py
+++ b/wetb/signal/fit/__init__.py
@@ -4,5 +4,6 @@ d = dir()
from wetb.signal.fit._linear_fit import *
from wetb.signal.fit._bin_fit import *
from wetb.signal.fit._fourier_fit import *
+from wetb.signal.fit._spline_fit import *
__all__ = sorted([m for m in set(dir()) - set(d)])
diff --git a/wetb/signal/fit/_bin_fit.py b/wetb/signal/fit/_bin_fit.py
index f1124ff6fb0919a22f2ef8a1e45ad14a597ce9eb..fcaf5879f9e0f89c99ded6cf8b69f32f48997e9e 100644
--- a/wetb/signal/fit/_bin_fit.py
+++ b/wetb/signal/fit/_bin_fit.py
@@ -103,14 +103,9 @@ def bin_fit(x,y, bins=10, kind='linear', bin_func=np.nanmean, bin_min_count=3, l
else:
raise NotImplementedError("Argument for handling upper observations, %s, not implemented"%upper)
- #Create mean function
- def fit(x):
- x = x[:].copy().astype(np.float)
- x[x<bin_x_fit[0]] = np.nan
- x[x>bin_x_fit[-1]] = np.nan
- return interp1d(bin_x_fit, bin_y_fit, kind)(x[:])
+
- return bin_x_fit, fit
+ return bin_x_fit, _interpolate_fit(bin_x_fit, bin_y_fit, kind)
def perpendicular_bin_fit(x, y, bins = 30, fit_func=None, bin_min_count=3, plt=None):
"""Fit a curve to the values, (x,y) using bins that are perpendicular to an initial fit
@@ -190,6 +185,15 @@ def perpendicular_bin_fit(x, y, bins = 30, fit_func=None, bin_min_count=3, plt=N
plt.plot(pbfx, pbfy, 'gray', label="perpendicular fit")
plt.legend()
#PlotData(None, bfx,bfy)
-
- return np.array(pc).T
+ bin_x_fit, bin_y_fit = np.array(pc).T
+ return bin_x_fit, _interpolate_fit(bin_x_fit, bin_y_fit, kind="linear")
+
+#Create mean function
+def _interpolate_fit(bin_x_fit, bin_y_fit, kind='linear'):
+ def fit(x):
+ x = x[:].copy().astype(np.float)
+ x[x<bin_x_fit[0]] = np.nan
+ x[x>bin_x_fit[-1]] = np.nan
+ return interp1d(bin_x_fit, bin_y_fit, kind)(x[:])
+ return fit
diff --git a/wetb/signal/fit/_spline_fit.py b/wetb/signal/fit/_spline_fit.py
new file mode 100644
index 0000000000000000000000000000000000000000..1fa074f9f359861797b42f3652fe8947fd938645
--- /dev/null
+++ b/wetb/signal/fit/_spline_fit.py
@@ -0,0 +1,74 @@
+import numpy as np
+def spline_fit(xp,yp):
+
+ def akima(x, y):
+ n = len(x)
+ var = np.zeros((n + 3))
+ z = np.zeros((n))
+ co = np.zeros((n, 4))
+ for i in range(n - 1):
+ var[i + 2] = (y[i + 1] - y[i]) / (x[i + 1] - x[i])
+ var[n + 1] = 2 * var[n] - var[n - 1]
+ var[n + 2] = 2 * var[n + 1] - var[n]
+ var[1] = 2 * var[2] - var[3]
+ var[0] = 2 * var[1] - var[2]
+
+ for i in range(n):
+ wi1 = abs(var[i + 3] - var[i + 2])
+ wi = abs(var[i + 1] - var[i])
+ if (wi1 + wi) == 0:
+ z[i] = (var[i + 2] + var[i + 1]) / 2
+ else:
+ z[i] = (wi1 * var[i + 1] + wi * var[i + 2]) / (wi1 + wi)
+
+ for i in range(n - 1):
+ dx = x[i + 1] - x[i]
+ a = (z[i + 1] - z[i]) * dx
+ b = y[i + 1] - y[i] - z[i] * dx
+ co[i, 0] = y[i]
+ co[i, 1] = z[i]
+ co[i, 2] = (3 * var[i + 2] - 2 * z[i] - z[i + 1]) / dx
+ co[i, 3] = (z[i] + z[i + 1] - 2 * var[i + 2]) / dx ** 2
+ co[n - 1, 0] = y[n - 1]
+ co[n - 1, 1] = z[n - 1]
+ co[n - 1, 2] = 0
+ co[n - 1, 3] = 0
+ return co
+
+ p_lst = [lambda x_, c=c, x0=x0: np.poly1d(c[::-1])(x_-x0) for c,x0 in zip(akima(xp,yp), xp)]
+
+ def spline(x):
+ y = np.empty_like(x)+np.nan
+ segment = np.searchsorted(xp,x, 'right')-1
+ for i in np.unique(segment):
+ m = segment==i
+ y[m] = p_lst[i](x[m])
+ return y
+# def coef2spline(x, xp, co):
+#
+# print (np.searchsorted(xp,x)-1)
+# x, y = [], []
+# for i, c in enumerate(co.tolist()[:-1]):
+# p = np.poly1d(c[::-1])
+# z = np.linspace(0, s[i + 1] - s[i ], 10, endpoint=i >= co.shape[0] - 2)
+# x.extend(s[i] + z)
+# y.extend(p(z))
+# return y
+#
+ return spline
+ #x, y, z = [coef2spline(curve_z_nd, akima(curve_z_nd, self.c2def[:, i])) for i in range(3)]
+ #return x, y, z
+
+if __name__=="__main__":
+ import matplotlib.pyplot as plt
+ x = np.random.randint(0,100,10)
+ t = np.arange(0,100,10)
+ plt.plot(t,x,'.',label='points')
+
+ t_ = np.arange(100)
+ spline = spline_fit(t,x)
+ print (np.abs(np.diff(np.diff(np.interp(t_, t,x)))).max())
+ print (np.abs(np.diff(np.diff(spline(t_)))).max())
+ plt.plot(t_, np.interp(t_, t,x))
+ plt.plot(t_, spline(t_),label='spline')
+ plt.show()
diff --git a/wetb/signal/fix/__init__.py b/wetb/signal/fix/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..d78316911251e8084e53f096c205f3f3620b9c22
--- /dev/null
+++ b/wetb/signal/fix/__init__.py
@@ -0,0 +1 @@
+from ._rotor_position import *
\ No newline at end of file
diff --git a/wetb/signal/fix/_rotor_position.py b/wetb/signal/fix/_rotor_position.py
new file mode 100644
index 0000000000000000000000000000000000000000..22aac76de54ba82a0cd4d21f501885e9a7364cb9
--- /dev/null
+++ b/wetb/signal/fix/_rotor_position.py
@@ -0,0 +1,113 @@
+'''
+Created on 30. mar. 2017
+
+@author: mmpe
+'''
+import numpy as np
+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
+ ----------
+ 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))
+ 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
+
+ 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)')
+
+ points = []
+ for j, i in enumerate(range(0,len(rp), N)):
+ #indexes for subsets for overlapping a and b polynomials
+ i1 = max(0,i-N2)
+ i2 = min(len(rp),i+N2)
+
+ #fit a polynomial
+ 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
+
+ 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)%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).
+ Optimal is defined as the value that minimizes the sum of squared differences
+ between differentiated rotor position and rotor speed
+
+ 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
+
+ """
+ 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)
+
+ best = 27
+ for step in [9,3,1]:
+
+ x_lst = np.arange(-2,3)*step + 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/interpolation.py b/wetb/signal/interpolation.py
new file mode 100644
index 0000000000000000000000000000000000000000..098a25b76f741eb215802ae0348807130b40e0b6
--- /dev/null
+++ b/wetb/signal/interpolation.py
@@ -0,0 +1,85 @@
+
+import numpy as np
+import numpy.ma as ma
+#from pylab import *
+def interpolate(x, xp, yp, max_xp_step=None, max_dydxp=None, cyclic_range=None, max_repeated=None):
+ """Interpolation similar to numpy.interp that handles nan and missing values
+
+ Parameters
+ ----------
+ x : array_like
+ The x-coordinates of the interpolated values.
+ xp : 1-D sequence of floats
+ The x-coordinates of the data points, must be increasing.
+ yp : 1-D sequence of floats
+ The y-coordinates of the data points, same length as xp.
+ max_xp_step : int, float or None, optional
+ Maximum xp-time step that is interpolated to x.\n
+ If time step > max_xp_step then NAN is returned for intermediate x values
+ If None, default, then this fix is not applied
+ max_dydxp : int, float, None, optional
+ Maximum absolute dydxp (slop of yp) that is interpolated to y.\n
+ If dydxp > max_dydxp then NAN is returned for intermediate y values
+ If None, default, then this fix is not applied
+ cyclick_range : int, float, None, optional
+ Range of posible values, e.g. 360 for degrees (both 0..360 and -180..180)
+ If None (default), data not interpreted as cyclic
+ max_repeated : int, float, None, optional
+ Maximum xp that yp are allowed to be repeated
+ if yp[i]==yp[i+1]==..==yp[i+j] and xp[i+j]-xp[i]>max_repeated_yp then
+ NAN is returned for xp[i]<x<=xp[i+j]
+
+
+ Returns
+ -------
+ y : {float, ndarray}
+ The interpolated values, same shape as x.
+
+ Examples
+ --------
+ >>> interpolate(x=[1,1.5,2,3], xp=[0.5,1.5,3], yp=[5,15,30])
+ [10,15,20,30]
+ >>> interpolate(x=[0, 1, 2, 7, 12], xp=[0, 2, 12], yp=[359, 0, 350], max_dydxp=45)
+ [359,nan,0,175,350]
+ """
+
+ xp = np.array(xp, dtype=np.float)
+ yp = np.array(yp, dtype=np.float)
+ assert xp.shape[0] == yp.shape[0], "xp and yp must have same length (%d!=%d)" % (xp.shape[0], yp.shape[0])
+ non_nan = ~(np.isnan(xp) & np.isnan(yp))
+ yp = yp[non_nan]
+ xp = xp[non_nan]
+ y = np.interp(x, xp, yp, np.nan, np.nan)
+
+
+ if cyclic_range is not None:
+ cr = cyclic_range
+ y2 = np.interp(x, xp, (yp + cr / 2) % cr - cr / 2, np.nan, np.nan) % cr
+ y = np.choose(np.r_[0, np.abs(np.diff(y)) > np.abs(np.diff(y2))], np.array([y, y2]))
+
+ if max_xp_step:
+ diff = np.diff(xp)
+ diff[np.isnan(diff)] = 0
+
+ indexes = (np.where(diff > max_xp_step)[0])
+ for index in indexes:
+ y[(x > xp[index]) & (x < xp[index + 1])] = np.nan
+ if max_dydxp:
+ if cyclic_range is None:
+ abs_dydxp = np.abs(np.diff(yp) / np.diff(xp))
+ else:
+ abs_dydxp = np.min([np.abs(np.diff((yp + cyclic_range / 2) % cyclic_range)) , np.abs(np.diff(yp % cyclic_range)) ], 0) / np.diff(xp)
+ abs_dydxp[np.isnan(abs_dydxp)] = 0
+
+ indexes = (np.where(abs_dydxp > max_dydxp)[0])
+ for index in indexes:
+ y[(x > xp[index]) & (x < xp[index + 1])] = np.nan
+ if max_repeated:
+ rep = np.r_[False, yp[1:] == yp[:-1], False]
+ tr = rep[1:] ^ rep[:-1]
+ itr = np.where(tr)[0]
+ for start, stop, l in zip (itr[::2] , itr[1::2], xp[itr[1::2]] - xp[itr[::2]]):
+ if l >= max_repeated:
+ y[(x > xp[start]) & (x <= xp[stop])] = np.nan
+ return y
+#print (interpolate(x=[0, 1, 2, 3, 4], xp=[0, 1, 2, 4], yp=[5, 5, 6, 6], max_repeated=1))
diff --git a/wetb/signal/nan_replace.py b/wetb/signal/nan_replace.py
new file mode 100644
index 0000000000000000000000000000000000000000..ffc691d0588eb7d22ce3894ea1dfe65b95ea48a1
--- /dev/null
+++ b/wetb/signal/nan_replace.py
@@ -0,0 +1,22 @@
+'''
+Created on 02/11/2015
+
+@author: MMPE
+'''
+
+import numpy as np
+from wetb.signal.filters import replacer
+def replace_by_mean(x):
+ return replacer.replace_by_mean(x, np.isnan(x))
+
+
+def replace_by_line(x):
+ return replacer.replace_by_line(x, np.isnan(x))
+
+def replace_by_polynomial(x, deg=3, no_base_points=12):
+ return replacer.replace_by_polynomial(x, np.isnan(x), deg, no_base_points)
+
+def max_no_nan(x):
+ return replacer.max_cont_mask_length(np.isnan(x))
+
+
diff --git a/wetb/signal/subset_mean.py b/wetb/signal/subset_mean.py
index 89ed1499cfa548a7c65403c444c6fb9a0f079dec..0de3ea896f8185e932998ddd85229b97523e2a9a 100644
--- a/wetb/signal/subset_mean.py
+++ b/wetb/signal/subset_mean.py
@@ -136,7 +136,53 @@ 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(values, rpm_dt=None, dmin=5, dmax=10, ):
+def revolution_trigger(rotor_position, sample_frq, rotor_speed, max_no_round_diff=1):
+ """Returns one index per revolution (minimum rotor position)
+
+ 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]
+
+ Returns
+ -------
+ nd_array : Array of indexes
+ """
+ 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]
+
+ 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
+ rp[np.r_[False, np.diff(rp)>thresshold]] = 180
+
+ upper_indexes = np.where((rp[:-1]>(360-thresshold))&(rp[1:]<(360-thresshold)))[0]
+ lower_indexes = np.where((rp[:-1]>thresshold)&(rp[1:]<thresshold))[0] +1
+
+ # 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
+
+
+def revolution_trigger_old(values, rpm_dt=None, dmin=5, dmax=10, ):
"""Return indexes where values are > max(values)-dmin and decreases more than dmax
If RPM and time step is provided, triggers steps < time of 1rpm is removed
diff --git a/wetb/signal/tests/test_differentiation.py b/wetb/signal/tests/test_differentiation.py
new file mode 100644
index 0000000000000000000000000000000000000000..87538cdc68ca0efa1a3ce733e57cc98f0ee1ab86
--- /dev/null
+++ b/wetb/signal/tests/test_differentiation.py
@@ -0,0 +1,20 @@
+'''
+Created on 29. mar. 2017
+
+@author: mmpe
+'''
+import unittest
+import numpy as np
+from wetb.signal.filters._differentiation import differentiation
+
+class Test(unittest.TestCase):
+
+
+ def testDifferentiation(self):
+ np.testing.assert_array_equal(differentiation([1,2,1,0,1,1]), [1,0,-1,0,.5,0])
+
+
+
+if __name__ == "__main__":
+ #import sys;sys.argv = ['', 'Test.testDifferentiation']
+ unittest.main()
\ No newline at end of file
diff --git a/wetb/signal/tests/test_fit.py b/wetb/signal/tests/test_fit.py
index ecdfd2cc69796b268d7351a376552d2d401a23ef..92c2c269a7759a9bfaeea492ff5e111c3d8a215d 100644
--- a/wetb/signal/tests/test_fit.py
+++ b/wetb/signal/tests/test_fit.py
@@ -10,6 +10,7 @@ import os
import unittest
from wetb.signal.fit import fourier_fit
from wetb.signal.error_measures import rms
+from wetb.signal.fit import spline_fit
tfp = os.path.join(os.path.dirname(__file__), 'test_files/')
class TestFit(unittest.TestCase):
@@ -88,11 +89,11 @@ class TestFit(unittest.TestCase):
x,y = ds('Wsp_metmast'), ds('Power')
if 0:
import matplotlib.pyplot as plt
- fx, fy = perpendicular_bin_fit(x,y,30,plt=plt)
+ fx, fit = perpendicular_bin_fit(x,y,30,plt=plt)
plt.show()
else:
- fx, fy = perpendicular_bin_fit(x,y,30)
+ fx, fit = perpendicular_bin_fit(x,y,30)
self.assertEqual(len(fx), 30)
@@ -143,6 +144,24 @@ class TestFit(unittest.TestCase):
# plt.plot(fourier_fit.F2x(np.fft.fft(y) / len(y)), label='fft')
# plt.legend()
# plt.show()
+
+ def test_spline(self):
+
+ x = np.random.randint(0,100,10)
+ t = np.arange(0,100,10)
+
+ t_ = np.arange(100)
+ spline = spline_fit(t,x)
+ acc_lin = np.diff(np.diff(np.interp(t_, t,x)))
+ acc_spline = np.diff(np.diff(spline(t_)))
+ self.assertLess(np.abs(acc_spline).max(), np.abs(acc_lin).max())
+ if 0:
+ import matplotlib.pyplot as plt
+ plt.plot(t,x,'.',label='points')
+ plt.plot(t_, spline(t_),label='spline')
+ plt.legend()
+ plt.show()
+
if __name__ == "__main__":
#import sys;sys.argv = ['', 'Test.testName']
unittest.main()
\ No newline at end of file
diff --git a/wetb/signal/tests/test_fix.py b/wetb/signal/tests/test_fix.py
new file mode 100644
index 0000000000000000000000000000000000000000..8663986285a6fdb9f3a37f487b751b619bd9cc83
--- /dev/null
+++ b/wetb/signal/tests/test_fix.py
@@ -0,0 +1,48 @@
+'''
+Created on 30. mar. 2017
+
+@author: mmpe
+'''
+import os
+import unittest
+from wetb import gtsdf
+from wetb.signal.fix._rotor_position import fix_rotor_position, find_fix_dt
+from wetb.signal.filters._differentiation import differentiation
+
+
+tfp = os.path.join(os.path.dirname(__file__), 'test_files/')
+import numpy as np
+class TestFix(unittest.TestCase):
+
+
+ def testRotorPositionFix(self):
+ ds = gtsdf.Dataset(tfp+'azi.hdf5')
+ sample_frq = 25
+
+ #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
+ err_sum = np.sum((rpm_pos - ds.Rot_cor)**2)
+
+ self.assertLess(err_sum,40)
+ if 0:
+ import matplotlib.pyplot as plt
+ 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()
+
+ def test_find_fix_dt(self):
+ ds = gtsdf.Dataset(tfp+'azi.hdf5')
+ sample_frq = 25
+
+ self.assertEqual(find_fix_dt(ds.azi, sample_frq, ds.Rot_cor), 4)
+
+
+if __name__ == "__main__":
+ #import sys;sys.argv = ['', 'Test.testRotorPositionFix']
+ unittest.main()
\ No newline at end of file
diff --git a/wetb/signal/tests/test_frq_filters.py b/wetb/signal/tests/test_frq_filters.py
new file mode 100644
index 0000000000000000000000000000000000000000..b4537de2c1f7986c67204db1046a22c1a5277247
--- /dev/null
+++ b/wetb/signal/tests/test_frq_filters.py
@@ -0,0 +1,94 @@
+'''
+Created on 27. mar. 2017
+
+@author: mmpe
+'''
+import unittest
+
+import numpy as np
+from wetb.signal.filters.frq_filters import sine_generator, low_pass, \
+ frequency_response, high_pass
+
+
+class Test(unittest.TestCase):
+
+
+ def test_sine_generator(self):
+ t,y = sine_generator(10,1,2)
+ self.assertEqual(np.diff(t).mean(),.1)
+ self.assertAlmostEqual(t.max(),1.9)
+ if 0:
+ import matplotlib.pyplot as plt
+ plt.plot(*sine_generator(10,1,2), label="1Hz sine")
+ plt.plot(*sine_generator(100,2,2), label="2Hz sine")
+ plt.legend()
+ plt.show()
+
+ def test_low_pass(self):
+ sf = 100 # sample frequency
+ t,y1 = sine_generator(sf,1,5) # 1Hz sine
+ t,y10 = sine_generator(sf,10,5) # 10Hz sine
+ sig = y1+y10
+ y_lp = low_pass(sig, sf, 1, order=1)
+
+ # 1 order:
+ # cut off frq: -3db
+ # Above cut off: 20db/decade
+ np.testing.assert_almost_equal(y_lp[100:400], y1[100:400]*.501, 1)
+ np.testing.assert_almost_equal((y_lp-y1*.501)[100:400], y10[100:400]*.01, 1)
+
+ if 0:
+ import matplotlib.pyplot as plt
+ plt.plot(t,sig, label='Input signal')
+ plt.plot(t, y1*0.501, label='1Hz sine (-3db)')
+ plt.plot(t,y_lp, label='Output signal')
+ plt.plot(t,y_lp - y1*0.501, label="Output - 1Hz sine(-3db)")
+ plt.plot(t, y10*0.01, label='10Hz sine (-20db)')
+ plt.plot(t, y_lp - y1*0.501 - y10*.01, label="(Output - 1Hz sine(-3db)) - 10Hz sine (-20db)")
+ plt.legend()
+ plt.show()
+
+ def test_frq_response(self):
+ w,h = frequency_response(100,10,'low',1)
+ self.assertAlmostEqual(np.interp(10, w, h), -3.01,2) # cut off frq -3.01 db
+ self.assertAlmostEqual(np.interp(100, w, h), -20,1) # -20 db per decade
+ if 0:
+ import matplotlib.pyplot as plt
+ frequency_response(100,10,'low',1, plt=plt)
+ frequency_response(100,10,'low',2, plt=plt)
+ frequency_response(100,10,'low',5, plt=plt)
+ frequency_response(100,10,'low',8, plt=plt)
+ frequency_response(100,10,'low',10, plt=plt)
+ frequency_response(100,10,'high',1, plt=plt)
+ frequency_response(100,10,'high',2, plt=plt)
+ plt.show()
+
+ def test_high_pass(self):
+ sf = 100 # sample frequency
+ t,y1 = sine_generator(sf,1,5) # 1Hz sine
+ t,y10 = sine_generator(sf,10,5) # 10Hz sine
+ sig = y1+y10
+ y_lp = high_pass(sig, sf, 10, order=1)
+
+ # 1 order:
+ # cut off frq: -3db
+ # Below cut off: 20db/decade
+ np.testing.assert_almost_equal(y_lp[100:400], y10[100:400]*.501, 1)
+ np.testing.assert_almost_equal((y_lp-y10*.501)[100:400], y1[100:400]*.01, 1)
+
+ if 0:
+ import matplotlib.pyplot as plt
+ plt.plot(t,sig, label='Input signal')
+ plt.plot(t, y10*0.501, label='10Hz sine (-3db)')
+ plt.plot(t,y_lp, label='Output signal')
+ plt.plot(t,y_lp - y10*0.501, label="Output - 10Hz sine(-3db)")
+ plt.plot(t, y1*0.01, label='1Hz sine (-20db)')
+ plt.plot(t, y_lp - y10*0.501 - y1*.01, label="(Output - 10Hz sine(-3db)) - 1Hz sine (-20db)")
+ plt.legend()
+ plt.show()
+
+
+
+if __name__ == "__main__":
+ #import sys;sys.argv = ['', 'Test.test_sine_generator']
+ unittest.main()
\ No newline at end of file
diff --git a/wetb/signal/tests/test_interpolation.py b/wetb/signal/tests/test_interpolation.py
new file mode 100644
index 0000000000000000000000000000000000000000..6fb71e8131a5d0f3765d25c154fd3910c2e410ca
--- /dev/null
+++ b/wetb/signal/tests/test_interpolation.py
@@ -0,0 +1,104 @@
+'''
+Created on 19/12/2014
+
+@author: MMPE
+'''
+import unittest
+
+from matplotlib.pyplot import *
+import numpy as np
+from wetb import signal
+from wetb import gtsdf
+import os
+from wetb.gtsdf.unix_time import from_unix
+import datetime
+
+
+class TestInterpolation(unittest.TestCase):
+
+ def test_interpolate1(self):
+ x = [1, 1.5, 2, 3]
+ xp = [0.5, 1.5, 3]
+ yp = [5, 15, 30]
+ np.testing.assert_array_equal(signal.interpolate(x, xp, yp), [10., 15., 20, 30.])
+ np.testing.assert_array_equal(signal.interpolate(x, xp, yp, 1), [10., 15., np.nan, 30.])
+
+
+
+
+ def test_interpolate2(self):
+ x = np.arange(0, 100, 5, dtype=np.float)
+ xp = np.arange(0, 100, 10, dtype=np.float)
+ xp = np.r_[xp[:3], xp[5:]]
+ yp = np.arange(10, dtype=np.float)
+ yp[7:8] = np.nan
+ yp = np.r_[yp[:3], yp[5:]]
+
+ #x = [ 0. 5. 10. 15. 20. 25. 30. 35. 40. 45. 50. 55. 60. 65. 70. 75. 80. 85. 90. 95.]
+ #xp = [0.0, 10.0, 20.0, 50.0, 60.0, 70.0, 80.0, 90.0]
+ #yp = [0.0, 1.0, 2.0, 5.0, 6.0, nan, 8.0, 9.0]
+
+ y = signal.interpolate(x, xp, yp)
+ np.testing.assert_array_equal(y[~np.isnan(y)], [0., 0.5, 1. , 1.5, 2. , 2.5, 3. , 3.5, 4. , 4.5, 5. , 5.5, 8. , 8.5, 9. ])
+
+ y = signal.interpolate(x, xp, yp, 10)
+ np.testing.assert_array_equal(y[~np.isnan(y)], [ 0. , 0.5, 1. , 1.5, 2. , 5. , 5.5, 8. , 8.5, 9. ])
+
+
+
+ def test_interpolate_max_dydxp(self):
+ x = np.arange(7)
+ xp = [0, 2, 4, 6]
+ yp = [358, 359, 0, 1]
+ y = signal.interpolate(x, xp, yp, max_dydxp=30)
+ np.testing.assert_array_equal(y, [ 358., 358.5, 359., np.nan, 0., 0.5, 1. ])
+
+ y = signal.interpolate(x, xp, yp, max_dydxp=180)
+ np.testing.assert_array_equal(y, [ 358., 358.5, 359., 179.5, 0., 0.5, 1. ])
+
+
+ def test_interpolate_max_dydxp_cyclic(self):
+ x = np.arange(7)
+ xp = [0, 2, 4, 6]
+ yp = [358, 359, 0, 1]
+
+ y = signal.interpolate(x, xp, [178, 179, -180, -179], max_dydxp=30, cyclic_range=360)
+ np.testing.assert_array_equal(y, [ 178. , 178.5, 179. , -0.5, -180. , -179.5, -179. ])
+
+ y = signal.interpolate(x, xp, yp, max_dydxp=30, cyclic_range=360)
+ np.testing.assert_array_equal(y, [ 358., 358.5, 359., 359.5, 0., 0.5, 1. ])
+
+ y = signal.interpolate(xp, x, [ 358., 358.5, 359., 359.5, 0., 0.5, 1. ], max_dydxp=30, cyclic_range=360)
+ np.testing.assert_array_equal(y, [ 358., 359., 0., 1. ])
+
+ y = signal.interpolate(x, xp, yp, max_dydxp=180)
+ np.testing.assert_array_equal(y, [ 358., 358.5, 359., 179.5, 0., 0.5, 1. ])
+
+ def test_interpolate_max_repeated(self):
+ x = np.arange(7)
+ xp = [0, 1, 2, 3, 4, 5, 6]
+ yp = [4, 5, 5, 5, 6, 6, 7]
+ y = signal.interpolate(x, xp, yp, max_repeated=2)
+ np.testing.assert_array_equal(y, [ 4, 5, np.nan, np.nan, 6, 6, 7])
+
+
+ x = np.arange(7)
+ xp = [0, 3, 4, 5, 6]
+ yp = [5, 5, 7, 6, 6]
+ y = signal.interpolate(x, xp, yp, max_repeated=2)
+ np.testing.assert_array_equal(y, [ 5, np.nan, np.nan, np.nan, 7, 6, 6])
+
+ xp = [0, 2, 3, 4, 6]
+ y = signal.interpolate(x, xp, yp, max_repeated=1)
+ np.testing.assert_array_equal(y, [ 5, np.nan, np.nan, 7, 6, np.nan, np.nan])
+
+
+ xp = [0, 1, 2, 3, 6]
+ y = signal.interpolate(x, xp, yp, max_repeated=2)
+ np.testing.assert_array_equal(y, [ 5, 5, 7, 6, np.nan, np.nan, np.nan])
+
+
+
+if __name__ == "__main__":
+ #import sys;sys.argv = ['', 'Test.testName']
+ unittest.main()
diff --git a/wetb/signal/tests/test_nan_replace.py b/wetb/signal/tests/test_nan_replace.py
new file mode 100644
index 0000000000000000000000000000000000000000..6211ce029a04ee604af3a3589f0df31ffe3838f5
--- /dev/null
+++ b/wetb/signal/tests/test_nan_replace.py
@@ -0,0 +1,66 @@
+'''
+Created on 02/11/2015
+
+@author: MMPE
+'''
+import unittest
+import numpy as np
+from wetb.signal.nan_replace import replace_by_mean, replace_by_line, \
+ replace_by_polynomial, max_no_nan
+from matplotlib.pyplot import plot, show
+class TestNan_replace(unittest.TestCase):
+
+
+ def test_nan_replace_by_mean(self):
+ a = np.array([1, 5, 6, 4, 5, np.nan, 3, 1, 5, 6, 4.])
+ np.testing.assert_array_equal(replace_by_mean(a), [1, 5, 6, 4, 5, 4, 3, 1, 5, 6, 4.])
+ a = np.array([1, 5, 6, 4, 5, np.nan, np.nan, 1, 5, 6, 4.])
+ np.testing.assert_array_equal(replace_by_mean(a), [1, 5, 6, 4, 5, 3, 3, 1, 5, 6, 4.])
+ a = np.array([np.nan, 5, 6, 4, 5, np.nan, 3, 1, 5, 6, np.nan])
+ np.testing.assert_array_equal(replace_by_mean(a), [5, 5, 6, 4, 5, 4, 3, 1, 5, 6, 6])
+
+
+ def test_nan_replace_by_line(self):
+ a = np.array([1, 5, 6, 4, 5, np.nan, 3, 1, 5, 6, 4.])
+ np.testing.assert_array_equal(replace_by_line(a), [1, 5, 6, 4, 5, 4, 3, 1, 5, 6, 4.])
+ a = np.array([1, 5, 6, 4, 5, np.nan, np.nan, 2, 5, 6, 4.])
+ np.testing.assert_array_equal(replace_by_line(a), [1, 5, 6, 4, 5, 4, 3, 2, 5, 6, 4.])
+ a = np.array([np.nan, 5, 6, 4, 5, np.nan, 3, 1, 5, 6, np.nan])
+ np.testing.assert_array_equal(replace_by_line(a), [5, 5, 6, 4, 5, 4, 3, 1, 5, 6, 6])
+
+
+ def test_nan_replace_by_polynomial(self):
+ a = np.array([np.nan, 5, 6, 4, 5, np.nan, 3, 1, 5, 6, np.nan])
+ np.testing.assert_array_equal(replace_by_polynomial(a), [5, 5, 6, 4, 5, 4, 3, 1, 5, 6, 6])
+ a = np.array([1, 5, 6, 4, 5, np.nan, np.nan, 1, 5, 6, 4.])
+ np.testing.assert_array_almost_equal(replace_by_polynomial(a, 3, 2), [1, 5, 6, 4, 5, 3.3, 1.2, 1, 5, 6, 4.])
+
+ if 0:
+ plot(a)
+ plot(replace_by_polynomial(a, 3, 2))
+ show()
+
+
+ def test_nan_replace_by_polynomial2(self):
+ a = np.r_[np.arange(10), np.repeat(np.nan, 100), 10 - np.arange(10)]
+ self.assertLessEqual(replace_by_polynomial(a, 3, 9).max(), np.nanmax(a))
+
+ if 0:
+ plot(a, '.r')
+ plot(replace_by_polynomial(a, 3, 9), 'g-')
+ show()
+
+
+ def test_max_no_nan(self):
+ a = np.array([1, 5, 6, 4, 5, np.nan, 3, 1, 5, np.nan, 4.])
+ self.assertEqual(max_no_nan(a), 1)
+ a = np.array([1, 5, 6, 4, 5, np.nan, np.nan, 1, 5, np.nan, 4.])
+ self.assertEqual(max_no_nan(a), 2)
+ a = np.array([np.nan, np.nan, np.nan, 4, 5, np.nan, np.nan, 1, 5, np.nan, 4.])
+ self.assertEqual(max_no_nan(a), 3)
+ a = np.array([1, 5, 6, 4, 5, np.nan, 4, 1, 5, np.nan, np.nan])
+ self.assertEqual(max_no_nan(a), 2)
+
+if __name__ == "__main__":
+ #import sys;sys.argv = ['', 'Test.test_nanreplace']
+ unittest.main()
diff --git a/wetb/signal/tests/test_subset_mean.py b/wetb/signal/tests/test_subset_mean.py
index 21f92445cba215d4e75357457e376fa85da1873a..23642ae7509bc80188499122807d5b22ddfce2fc 100644
--- a/wetb/signal/tests/test_subset_mean.py
+++ b/wetb/signal/tests/test_subset_mean.py
@@ -9,7 +9,7 @@ import unittest
import numpy as np
from wetb import gtsdf
from wetb.signal.subset_mean import time_trigger, subset_mean, \
- non_nan_index_trigger, revolution_trigger
+ non_nan_index_trigger, revolution_trigger_old, revolution_trigger
from wetb.utils.geometry import rpm2rads
@@ -68,10 +68,10 @@ class TestSubsetMean(unittest.TestCase):
ds = gtsdf.Dataset(tfp+'azi.hdf5')
azi, rpm, time = [ds(x)[8403:8803] for x in ['azi','Rot_cor','Time']]
- trigger = revolution_trigger(azi)
+ trigger = revolution_trigger_old(azi)
np.testing.assert_array_equal(trigger, [ 17, 128, 241, 354])
azi[64] = 358
- trigger = revolution_trigger(azi, (ds('Rot_cor'), np.diff(time).mean()))
+ trigger = revolution_trigger_old(azi, (ds('Rot_cor'), np.diff(time).mean()))
# import matplotlib.pyplot as plt
# t = np.arange(len(azi))
@@ -82,5 +82,31 @@ class TestSubsetMean(unittest.TestCase):
np.testing.assert_array_equal(trigger, [ (128,241),(241,354)])
+ def test_revolution_trigger(self):
+ rotor_position = np.arange(0.,360*10,4)
+ rotor_position += np.random.random(len(rotor_position))
+ rotor_position = rotor_position % 360
+ if 0:
+ x1 = np.random.randint(0, len(rotor_position),10)
+ print (list(x1))
+ x2 = np.random.randint(0, len(rotor_position),10)
+ print (list(x2))
+ else:
+ x1 = [447, 854, 595, 804, 847, 488, 412, 199, 675, 766]
+ x2 = [92, 647, 821, 422, 33, 159, 369, 99, 157, 464]
+
+ rotor_position[x1] += 360
+ rotor_position[x2] -= 360
+ rotor_position[90] = 180
+
+ 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()
+
+
if __name__ == "__main__":
unittest.main()
diff --git a/wetb/wind/shear.py b/wetb/wind/shear.py
index dd6122470e3a5a460660b7673371fb9042a147dd..dd3281b0ff5139253d76487a7be6893ad538e6f2 100644
--- a/wetb/wind/shear.py
+++ b/wetb/wind/shear.py
@@ -99,7 +99,7 @@ def fit_power_shear_ref(z_u_lst, z_ref, plt=None):
"""
def shear_error(x, z_u_lst, z_ref):
alpha, u_ref = x
- return np.sum([(u - u_ref * (z / z_ref) ** alpha) ** 2 for z, u in z_u_lst])
+ return np.nansum([(u - u_ref * (z / z_ref) ** alpha) ** 2 for z, u in z_u_lst])
z_u_lst = [(z, np.mean(u)) for z, u in z_u_lst]
alpha, u_ref = fmin(shear_error, (.1, 10), (z_u_lst, z_ref), disp=False)
if plt:
@@ -107,6 +107,9 @@ def fit_power_shear_ref(z_u_lst, z_ref, plt=None):
plt.plot(u, z, '.')
z = np.linspace(min(z), max(z), 100)
plt.plot(power_shear(alpha, z_ref, u_ref)(z), z)
+ plt.margins(.1)
+ if alpha==.1 and u_ref==10: # Initial conditions
+ return np.nan, np.nan
return alpha, u_ref
diff --git a/wetb/wind/tests/test_Shear.py b/wetb/wind/tests/test_Shear.py
index 38f1d0db03c0148b86911304d103f607aa43108e..44a369a2cbac3cb665138272089a69f242935d7e 100644
--- a/wetb/wind/tests/test_Shear.py
+++ b/wetb/wind/tests/test_Shear.py
@@ -49,6 +49,24 @@ class TestShear(unittest.TestCase):
14: WSP gl. coo.,Vz 21
"""
+ def test_power_shear_fit(self):
+ z = [30,50,70]
+ a,u_ref = .3,10
+ u = power_shear(a,50,u_ref)(z)
+ np.testing.assert_array_almost_equal(fit_power_shear_ref(zip(z,u), 50), [a,u_ref],3)
+ if 0:
+ import matplotlib.pyplot as plt
+ fit_power_shear_ref(zip(z,u), 50, plt)
+ plt.show()
+
+ def test_power_shear_fit_nan(self):
+ z = [30,50,70]
+ a,u_ref = .3,10
+ u = power_shear(a,50,u_ref)(z)
+ u[2] = np.nan
+ np.testing.assert_array_almost_equal(fit_power_shear_ref(zip(z,u), 50), [a,u_ref],3)
+ u[:] = np.nan
+ np.testing.assert_array_almost_equal(fit_power_shear_ref(zip(z,u), 50), [np.nan, np.nan],3)
def test_power_shear(self):
if all: