diff --git a/README.md b/README.md
index 1900083802a855519b34d211dbd7a09836ce87a5..2365c88119f9fa301de7f6691503695292ed49fe 100644
--- a/README.md
+++ b/README.md
@@ -25,6 +25,7 @@ General Time Series Data Format, a binary hdf5 data format for storing time seri
### [fatigue_tools](wetb/fatigue_tools)
- [fatigue](wetb/fatigue_tools/fatigue.py): Rainflow counting, cycle matrix and equvivalent loads
+- [bearing_damage](wetb/fatigue_tools/bearing_damage.py): Calculate a comparable measure of bearing damage
### [wind](wetb/wind)
- [shear](wetb/wind/shear.py): Calculate and fit wind shear
diff --git a/wetb/fatigue_tools/fatigue.py b/wetb/fatigue_tools/fatigue.py
index 0c7425f4d790b559a992eacaacf0cfb4febed2d2..a411a1674a222fc7a6066970b225ec893b5f5eb7 100644
--- a/wetb/fatigue_tools/fatigue.py
+++ b/wetb/fatigue_tools/fatigue.py
@@ -2,199 +2,23 @@
Created on 04/03/2013
@author: mmpe
-'rainflow_windap' and 'rainflow_astm' are two different methods to for rainflow counting
+
'eq_load' calculate equivalent loads using one of the two rain flow counting methods
'cycle_matrix' calculates a matrix of cycles (binned on amplitude and mean value)
'eq_load_and_cycles' is used to calculate eq_loads of multiple time series (e.g. life time equivalent load)
+The methods uses the rainflow counting routines (See documentation in top of methods):
+- 'rainflow_windap': (Described in "Recommended Practices for Wind Turbine Testing - 3. Fatigue Loads",
+ 2. edition 1990, Appendix A)
+or
+- 'rainflow_astm' (based on the c-implementation by Adam Nieslony found at the MATLAB Central File Exchange
+ http://www.mathworks.com/matlabcentral/fileexchange/3026)
'''
-
-#try:
-# """
-# The cython_import function compiles modules using cython.
-# It is found at: https://github.com/madsmpedersen/MMPE/blob/master/cython_compile/cython_compile.py
-# """
-# from mmpe.cython_compile.cython_compile import cython_import
-#except ImportError:
-# cython_import = __import__
import numpy as np
+from wetb.fatigue_tools.rainflowcounting import rainflowcount
-
-def rfc_hist(sig_rf, nrbins=46):
- """
- Histogram of rainflow counted cycles
- ====================================
-
- hist, bin_edges, bin_avg = rfc_hist(sig, nrbins=46)
-
- Divide the rainflow counted cycles of a signal into equally spaced bins.
-
- Created on Wed Feb 16 16:53:18 2011
- @author: David Verelst
- Modified 10.10.2011 by Mads M Pedersen to elimintate __copy__ and __eq__
-
- Parameters
- ----------
- sig_rf : array-like
- As output by rfc_astm or rainflow
-
- nrbins : int, optional
- Divide the rainflow counted amplitudes in a number of equally spaced
- bins.
-
- Returns
- -------
- hist : array-like
- Counted rainflow cycles per bin, has nrbins elements
-
- bin_edges : array-like
- Edges of the bins, has nrbins+1 elements.
-
- bin_avg : array-like
- Average rainflow cycle amplitude per bin, has nrbins elements.
- """
-
- rf_half = sig_rf
-
- # the Matlab approach is to divide into 46 bins
- bin_edges = np.linspace(0, 1, num=nrbins + 1) * rf_half.max()
- hist = np.histogram(rf_half, bins=bin_edges)[0]
- # calculate the average per bin
- hist_sum = np.histogram(rf_half, weights=rf_half, bins=bin_edges)[0]
- # replace zeros with one, to avoid 0/0
- hist_ = hist.copy()
- hist_[(hist == 0).nonzero()] = 1.0
- # since the sum is also 0, the avg remains zero for those whos hist is zero
- bin_avg = hist_sum / hist_
-
- return hist, bin_edges, bin_avg
-
-
-def check_signal(signal):
- # check input data validity
- if not type(signal).__name__ == 'ndarray':
- raise TypeError('signal must be ndarray, not: ' + type(signal).__name__)
-
- elif len(signal.shape) not in (1, 2):
- raise TypeError('signal must be 1D or 2D, not: ' + str(len(signal.shape)))
-
- if len(signal.shape) == 2:
- if signal.shape[1] > 1:
- raise TypeError('signal must have one column only, not: ' + str(signal.shape[1]))
- if np.min(signal) == np.max(signal):
- raise TypeError("Signal contains no variation")
-
-
-def rainflow_windap(signal, levels=255., thresshold=(255 / 50)):
- """Windap equivalent rainflow counting
-
-
- Calculate the amplitude and mean values of half cycles in signal
-
- This algorithms used by this routine is implemented directly as described in
- "Recommended Practices for Wind Turbine Testing - 3. Fatigue Loads", 2. edition 1990, Appendix A
-
- Parameters
- ----------
- Signal : array-like
- The raw signal
-
- levels : int, optional
- The signal is discretize into this number of levels.
- 255 is equivalent to the implementation in Windap
-
- thresshold : int, optional
- Cycles smaller than this thresshold are ignored
- 255/50 is equivalent to the implementation in Windap
-
- Returns
- -------
- ampl : array-like
- Peak to peak amplitudes of the half cycles
-
- mean : array-like
- Mean values of the half cycles
-
-
- Example
- -------
- >>> signal = np.array([-2.0, 0.0, 1.0, 0.0, -3.0, 0.0, 5.0, 0.0, -1.0, 0.0, 3.0, 0.0, -4.0, 0.0, 4.0, 0.0, -2.0])
- >>> ampl, mean = rainflow_windap(signal)
- """
- check_signal(signal)
- #type <double> is required by <find_extreme> and <rainflow>
- signal = signal.astype(np.double)
-
- offset = np.nanmin(signal)
- signal -= offset
- if np.nanmax(signal) > 0:
- gain = np.nanmax(signal) / levels
- signal = signal / gain
- signal = np.round(signal).astype(np.int)
-
- from wetb.fatigue_tools.peak_trough import peak_trough
-
-
- #Convert to list of local minima/maxima where difference > thresshold
- sig_ext = peak_trough(signal, thresshold)
-
- from wetb.fatigue_tools.pair_range import pair_range_amplitude_mean
-
- #rainflow count
- ampl_mean = pair_range_amplitude_mean(sig_ext)
-
- ampl_mean = np.array(ampl_mean)
- ampl_mean = np.round(ampl_mean / thresshold) * gain * thresshold
- ampl_mean[:, 1] += offset
- return ampl_mean.T
-
-
-
-def rainflow_astm(signal):
- """Matlab equivalent rainflow counting
-
- Calculate the amplitude and mean values of half cycles in signal
-
- This implemementation is based on the c-implementation by Adam Nieslony found at
- the MATLAB Central File Exchange http://www.mathworks.com/matlabcentral/fileexchange/3026
-
- Parameters
- ----------
- Signal : array-like
- The raw signal
-
- Returns
- -------
- ampl : array-like
- peak to peak amplitudes of the half cycles (note that the matlab implementation
- uses peak amplitude instead of peak to peak)
-
- mean : array-like
- Mean values of the half cycles
-
-
- Examples
- --------
- >>> signal = np.array([-2.0, 0.0, 1.0, 0.0, -3.0, 0.0, 5.0, 0.0, -1.0, 0.0, 3.0, 0.0, -4.0, 0.0, 4.0, 0.0, -2.0])
- >>> ampl, mean = rainflow_astm(signal)
- """
- check_signal(signal)
-
- # type <double> is reuqired by <find_extreme> and <rainflow>
- signal = signal.astype(np.double)
-
- # Import find extremes and rainflow.
- # If possible the module is compiled using cython otherwise the python implementation is used
- #cython_import('fatigue_tools.rainflowcount_astm')
- from wetb.fatigue_tools.rainflowcount_astm import find_extremes, rainflowcount
-
- # Remove points which is not local minimum/maximum
- sig_ext = find_extremes(signal)
-
- # rainflow count
- ampl_mean = np.array(rainflowcount(sig_ext))
-
- return np.array(ampl_mean).T
+rainflow_windap = rainflowcount.rainflow_windap
+rainflow_astm = rainflowcount.rainflow_astm
def eq_load(signals, no_bins=46, m=[3, 4, 6, 8, 10, 12], neq=1, rainflow_func=rainflow_windap):
@@ -239,14 +63,6 @@ def eq_load(signals, no_bins=46, m=[3, 4, 6, 8, 10, 12], neq=1, rainflow_func=ra
except TypeError:
return [[np.nan] * len(np.atleast_1d(m))] * len(np.atleast_1d(neq))
-# ampl, _ = rainflow_func(signals)
-# if ampl is None:
-# return []
-# hist_data, x, bin_avg = rfc_hist(ampl, no_bins)
-#
-# m = np.atleast_1d(m)
-#
-# return np.array([np.power(np.sum(0.5 * hist_data * np.power(bin_avg, m[i])) / neq, 1. / m[i]) for i in range(len(m))])
def eq_load_and_cycles(signals, no_bins=46, m=[3, 4, 6, 8, 10, 12], neq=[10 ** 6, 10 ** 7, 10 ** 8], rainflow_func=rainflow_windap):
@@ -354,3 +170,18 @@ def cycle_matrix(signals, ampl_bins=10, mean_bins=10, rainflow_func=rainflow_win
return cycles, ampl_bin_mean, ampl_edges, mean_bin_mean, mean_edges
+if __name__ == "__main__":
+ signal1 = np.array([-2.0, 0.0, 1.0, 0.0, -3.0, 0.0, 5.0, 0.0, -1.0, 0.0, 3.0, 0.0, -4.0, 0.0, 4.0, 0.0, -2.0])
+ signal2 = signal1 * 1.1
+
+ # equivalent load for default wöhler slopes
+ print (eq_load(signal1, no_bins=50, neq=17, rainflow_func=rainflow_windap))
+ print (eq_load(signal1, no_bins=50, neq=17, rainflow_func=rainflow_astm))
+
+ # Cycle matrix with 4 amplitude bins and 4 mean value bins
+ print (cycle_matrix(signal1, 4, 4, rainflow_func=rainflow_windap))
+ print (cycle_matrix(signal1, 4, 4, rainflow_func=rainflow_astm))
+
+ # Cycle matrix where signal1 and signal2 contributes with 50% each
+ print (cycle_matrix([(.5, signal1), (.5, signal2)], 4, 8, rainflow_func=rainflow_astm))
+
diff --git a/wetb/fatigue_tools/rainflowcounting/__init__.py b/wetb/fatigue_tools/rainflowcounting/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/wetb/fatigue_tools/rainflowcounting/compile.py b/wetb/fatigue_tools/rainflowcounting/compile.py
new file mode 100644
index 0000000000000000000000000000000000000000..26663108f33745e6c1f5dbad1829cc0abb992ffd
--- /dev/null
+++ b/wetb/fatigue_tools/rainflowcounting/compile.py
@@ -0,0 +1,7 @@
+import sys
+sys.path.append("../../../../MMPE/")
+from mmpe.cython_compile.cython_compile import cython_import
+
+cython_import('pair_range')
+cython_import('peak_trough')
+cython_import('rainflowcount_astm')
diff --git a/wetb/fatigue_tools/pair_range.py b/wetb/fatigue_tools/rainflowcounting/pair_range.py
similarity index 100%
rename from wetb/fatigue_tools/pair_range.py
rename to wetb/fatigue_tools/rainflowcounting/pair_range.py
diff --git a/wetb/fatigue_tools/pair_range.pyd b/wetb/fatigue_tools/rainflowcounting/pair_range.pyd
similarity index 71%
rename from wetb/fatigue_tools/pair_range.pyd
rename to wetb/fatigue_tools/rainflowcounting/pair_range.pyd
index e3e28cb4e3bb7ab0fc866c3ab7d1ec242458abdb..30f8baed7a441394607bf0e5759a496a08e83a6d 100644
Binary files a/wetb/fatigue_tools/pair_range.pyd and b/wetb/fatigue_tools/rainflowcounting/pair_range.pyd differ
diff --git a/wetb/fatigue_tools/peak_trough.py b/wetb/fatigue_tools/rainflowcounting/peak_trough.py
similarity index 100%
rename from wetb/fatigue_tools/peak_trough.py
rename to wetb/fatigue_tools/rainflowcounting/peak_trough.py
diff --git a/wetb/fatigue_tools/peak_trough.pyd b/wetb/fatigue_tools/rainflowcounting/peak_trough.pyd
similarity index 99%
rename from wetb/fatigue_tools/peak_trough.pyd
rename to wetb/fatigue_tools/rainflowcounting/peak_trough.pyd
index db4087cd89878143c3de4f75a60b5026f5dceb1f..250f85f0e8de0b2c2fe387e0175e92e9c8415cf7 100644
Binary files a/wetb/fatigue_tools/peak_trough.pyd and b/wetb/fatigue_tools/rainflowcounting/peak_trough.pyd differ
diff --git a/wetb/fatigue_tools/rainflowcounting/rainflowcount.py b/wetb/fatigue_tools/rainflowcounting/rainflowcount.py
new file mode 100644
index 0000000000000000000000000000000000000000..f05a1a0a1f38f4fe17c3a436a8be7e548e0e133a
--- /dev/null
+++ b/wetb/fatigue_tools/rainflowcounting/rainflowcount.py
@@ -0,0 +1,126 @@
+import numpy as np
+def check_signal(signal):
+ # check input data validity
+ if not type(signal).__name__ == 'ndarray':
+ raise TypeError('signal must be ndarray, not: ' + type(signal).__name__)
+
+ elif len(signal.shape) not in (1, 2):
+ raise TypeError('signal must be 1D or 2D, not: ' + str(len(signal.shape)))
+
+ if len(signal.shape) == 2:
+ if signal.shape[1] > 1:
+ raise TypeError('signal must have one column only, not: ' + str(signal.shape[1]))
+ if np.min(signal) == np.max(signal):
+ raise TypeError("Signal contains no variation")
+
+
+def rainflow_windap(signal, levels=255., thresshold=(255 / 50)):
+ """Windap equivalent rainflow counting
+
+
+ Calculate the amplitude and mean values of half cycles in signal
+
+ This algorithms used by this routine is implemented directly as described in
+ "Recommended Practices for Wind Turbine Testing - 3. Fatigue Loads", 2. edition 1990, Appendix A
+
+ Parameters
+ ----------
+ Signal : array-like
+ The raw signal
+
+ levels : int, optional
+ The signal is discretize into this number of levels.
+ 255 is equivalent to the implementation in Windap
+
+ thresshold : int, optional
+ Cycles smaller than this thresshold are ignored
+ 255/50 is equivalent to the implementation in Windap
+
+ Returns
+ -------
+ ampl : array-like
+ Peak to peak amplitudes of the half cycles
+
+ mean : array-like
+ Mean values of the half cycles
+
+
+ Example
+ -------
+ >>> signal = np.array([-2.0, 0.0, 1.0, 0.0, -3.0, 0.0, 5.0, 0.0, -1.0, 0.0, 3.0, 0.0, -4.0, 0.0, 4.0, 0.0, -2.0])
+ >>> ampl, mean = rainflow_windap(signal)
+ """
+ check_signal(signal)
+ #type <double> is required by <find_extreme> and <rainflow>
+ signal = signal.astype(np.double)
+
+ offset = np.nanmin(signal)
+ signal -= offset
+ if np.nanmax(signal) > 0:
+ gain = np.nanmax(signal) / levels
+ signal = signal / gain
+ signal = np.round(signal).astype(np.int)
+
+ from wetb.fatigue_tools.rainflowcounting.peak_trough import peak_trough
+
+
+ #Convert to list of local minima/maxima where difference > thresshold
+ sig_ext = peak_trough(signal, thresshold)
+
+ from wetb.fatigue_tools.rainflowcounting.pair_range import pair_range_amplitude_mean
+
+ #rainflow count
+ ampl_mean = pair_range_amplitude_mean(sig_ext)
+
+ ampl_mean = np.array(ampl_mean)
+ ampl_mean = np.round(ampl_mean / thresshold) * gain * thresshold
+ ampl_mean[:, 1] += offset
+ return ampl_mean.T
+
+
+
+def rainflow_astm(signal):
+ """Matlab equivalent rainflow counting
+
+ Calculate the amplitude and mean values of half cycles in signal
+
+ This implemementation is based on the c-implementation by Adam Nieslony found at
+ the MATLAB Central File Exchange http://www.mathworks.com/matlabcentral/fileexchange/3026
+
+ Parameters
+ ----------
+ Signal : array-like
+ The raw signal
+
+ Returns
+ -------
+ ampl : array-like
+ peak to peak amplitudes of the half cycles (note that the matlab implementation
+ uses peak amplitude instead of peak to peak)
+
+ mean : array-like
+ Mean values of the half cycles
+
+
+ Examples
+ --------
+ >>> signal = np.array([-2.0, 0.0, 1.0, 0.0, -3.0, 0.0, 5.0, 0.0, -1.0, 0.0, 3.0, 0.0, -4.0, 0.0, 4.0, 0.0, -2.0])
+ >>> ampl, mean = rainflow_astm(signal)
+ """
+ check_signal(signal)
+
+ # type <double> is reuqired by <find_extreme> and <rainflow>
+ signal = signal.astype(np.double)
+
+ # Import find extremes and rainflow.
+ # If possible the module is compiled using cython otherwise the python implementation is used
+ #cython_import('fatigue_tools.rainflowcount_astm')
+ from wetb.fatigue_tools.rainflowcounting.rainflowcount_astm import find_extremes, rainflowcount
+
+ # Remove points which is not local minimum/maximum
+ sig_ext = find_extremes(signal)
+
+ # rainflow count
+ ampl_mean = np.array(rainflowcount(sig_ext))
+
+ return np.array(ampl_mean).T
\ No newline at end of file
diff --git a/wetb/fatigue_tools/rainflowcount_astm.py b/wetb/fatigue_tools/rainflowcounting/rainflowcount_astm.py
similarity index 100%
rename from wetb/fatigue_tools/rainflowcount_astm.py
rename to wetb/fatigue_tools/rainflowcounting/rainflowcount_astm.py
diff --git a/wetb/fatigue_tools/rainflowcount_astm.pyd b/wetb/fatigue_tools/rainflowcounting/rainflowcount_astm.pyd
similarity index 99%
rename from wetb/fatigue_tools/rainflowcount_astm.pyd
rename to wetb/fatigue_tools/rainflowcounting/rainflowcount_astm.pyd
index 6ef6cd427dfb019332137805ebe4db80d90dcbc8..e985824ec83011c3b893870352cdef19cba7a582 100644
Binary files a/wetb/fatigue_tools/rainflowcount_astm.pyd and b/wetb/fatigue_tools/rainflowcounting/rainflowcount_astm.pyd differ
diff --git a/wetb/fatigue_tools/rainflowcounting/rfc_hist.py b/wetb/fatigue_tools/rainflowcounting/rfc_hist.py
new file mode 100644
index 0000000000000000000000000000000000000000..d4c7fb0a796e284d760239f764ffbe9e994441fc
--- /dev/null
+++ b/wetb/fatigue_tools/rainflowcounting/rfc_hist.py
@@ -0,0 +1,49 @@
+import numpy as np
+def rfc_hist(sig_rf, nrbins=46):
+ """
+ Histogram of rainflow counted cycles
+ ====================================
+
+ hist, bin_edges, bin_avg = rfc_hist(sig, nrbins=46)
+
+ Divide the rainflow counted cycles of a signal into equally spaced bins.
+
+ Created on Wed Feb 16 16:53:18 2011
+ @author: David Verelst
+ Modified 10.10.2011 by Mads M Pedersen to elimintate __copy__ and __eq__
+
+ Parameters
+ ----------
+ sig_rf : array-like
+ As output by rfc_astm or rainflow
+
+ nrbins : int, optional
+ Divide the rainflow counted amplitudes in a number of equally spaced
+ bins.
+
+ Returns
+ -------
+ hist : array-like
+ Counted rainflow cycles per bin, has nrbins elements
+
+ bin_edges : array-like
+ Edges of the bins, has nrbins+1 elements.
+
+ bin_avg : array-like
+ Average rainflow cycle amplitude per bin, has nrbins elements.
+ """
+
+ rf_half = sig_rf
+
+ # the Matlab approach is to divide into 46 bins
+ bin_edges = np.linspace(0, 1, num=nrbins + 1) * rf_half.max()
+ hist = np.histogram(rf_half, bins=bin_edges)[0]
+ # calculate the average per bin
+ hist_sum = np.histogram(rf_half, weights=rf_half, bins=bin_edges)[0]
+ # replace zeros with one, to avoid 0/0
+ hist_ = hist.copy()
+ hist_[(hist == 0).nonzero()] = 1.0
+ # since the sum is also 0, the avg remains zero for those whos hist is zero
+ bin_avg = hist_sum / hist_
+
+ return hist, bin_edges, bin_avg
\ No newline at end of file