Skip to content
Snippets Groups Projects
Commit 00c83d7a authored by Mads M. Pedersen's avatar Mads M. Pedersen
Browse files

test similar to matlab test at https://se.mathworks.com/help/signal/ref/rainflow.html

parent 694dc264
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
wetb/fatigue_tools/tests/test_fatigue.py
+ 51
− 32
View file @ 00c83d7a
......@@ -21,6 +21,7 @@ import os
testfilepath = os.path.join(os.path.dirname(__file__), 'test_files/') # test file path
class TestFatigueTools(unittest.TestCase):
def test_leq_1hz(self):
......@@ -31,26 +32,26 @@ class TestFatigueTools(unittest.TestCase):
m = 1
point_per_deg = 100
for amplitude in [1,2,3]:
for amplitude in [1, 2, 3]:
peak2peak = amplitude * 2
# sine signal with 10 periods (20 peaks)
nr_periods = 10
time = np.linspace(0, nr_periods*2*np.pi, point_per_deg*180)
time = np.linspace(0, nr_periods * 2 * np.pi, point_per_deg * 180)
neq = time[-1]
# mean value of the signal shouldn't matter
signal = amplitude * np.sin(time) + 5
r_eq_1hz = eq_load(signal, no_bins=1, m=m, neq=neq)[0]
r_eq_1hz_expected = ((2*nr_periods*amplitude**m)/neq)**(1/m)
r_eq_1hz_expected = ((2 * nr_periods * amplitude**m) / neq)**(1 / m)
np.testing.assert_allclose(r_eq_1hz, r_eq_1hz_expected)
# sine signal with 20 periods (40 peaks)
nr_periods = 20
time = np.linspace(0, nr_periods*2*np.pi, point_per_deg*180)
time = np.linspace(0, nr_periods * 2 * np.pi, point_per_deg * 180)
neq = time[-1]
# mean value of the signal shouldn't matter
signal = amplitude * np.sin(time) + 9
r_eq_1hz2 = eq_load(signal, no_bins=1, m=m, neq=neq)[0]
r_eq_1hz_expected2 = ((2*nr_periods*amplitude**m)/neq)**(1/m)
r_eq_1hz_expected2 = ((2 * nr_periods * amplitude**m) / neq)**(1 / m)
np.testing.assert_allclose(r_eq_1hz2, r_eq_1hz_expected2)
# 1hz equivalent should be independent of the length of the signal
......@@ -66,46 +67,44 @@ class TestFatigueTools(unittest.TestCase):
point_per_deg = 100
nr_periods = 10
time = np.linspace(0, nr_periods*2*np.pi, point_per_deg*180)
time = np.linspace(0, nr_periods * 2 * np.pi, point_per_deg * 180)
signal = (amplitude*np.sin(time)) + 5 + (amplitude*0.2*np.cos(5*time))
signal = (amplitude * np.sin(time)) + 5 + (amplitude * 0.2 * np.cos(5 * time))
cycles, ampl_bin_mean, ampl_edges, mean_bin_mean, mean_edges = \
cycle_matrix(signal, ampl_bins=10, mean_bins=5)
cycles.sum()
def test_astm1(self):
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)
np.testing.assert_array_equal(np.histogram2d(ampl, mean, [6, 4])[0], np.array([[ 0., 1., 0., 0.],
[ 1., 0., 0., 2.],
[ 0., 0., 0., 0.],
[ 0., 0., 0., 1.],
[ 0., 0., 0., 0.],
[ 0., 0., 1., 2.]]))
np.testing.assert_array_equal(np.histogram2d(ampl, mean, [6, 4])[0], np.array([[0., 1., 0., 0.],
[1., 0., 0., 2.],
[0., 0., 0., 0.],
[0., 0., 0., 1.],
[0., 0., 0., 0.],
[0., 0., 1., 2.]]))
def test_windap1(self):
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, 18, 2)
np.testing.assert_array_equal(np.histogram2d(ampl, mean, [6, 4])[0], np.array([[ 0., 0., 1., 0.],
[ 1., 0., 0., 2.],
[ 0., 0., 0., 0.],
[ 0., 0., 0., 1.],
[ 0., 0., 0., 0.],
[ 0., 0., 2., 1.]]))
np.testing.assert_array_equal(np.histogram2d(ampl, mean, [6, 4])[0], np.array([[0., 0., 1., 0.],
[1., 0., 0., 2.],
[0., 0., 0., 0.],
[0., 0., 0., 1.],
[0., 0., 0., 0.],
[0., 0., 2., 1.]]))
def test_windap2(self):
data = Hawc2io.ReadHawc2(testfilepath + "test").ReadBinary([2]).flatten()
np.testing.assert_allclose(eq_load(data, neq=61), np.array([[1.356, 1.758, 2.370, 2.784, 3.077, 3.296]]), 0.01)
def test_astm2(self):
data = Hawc2io.ReadHawc2(testfilepath + "test").ReadBinary([2]).flatten()
np.testing.assert_allclose(eq_load(data, neq=61, rainflow_func=rainflow_astm), np.array([[1.356, 1.758, 2.370, 2.784, 3.077, 3.296]]), 0.01)
np.testing.assert_allclose(eq_load(data, neq=61, rainflow_func=rainflow_astm),
np.array([[1.356, 1.758, 2.370, 2.784, 3.077, 3.296]]), 0.01)
# def test_windap3(self):
......@@ -118,20 +117,40 @@ class TestFatigueTools(unittest.TestCase):
# [ 0., 0., 0., 0.],
# [ 0., 1., 2., 0.]]) / 2)
def test_astm3(self):
data = Hawc2io.ReadHawc2(testfilepath + "test").ReadBinary([2]).flatten()
np.testing.assert_allclose(cycle_matrix(data, 4, 4, rainflow_func=rainflow_astm)[0], np.array([[ 24., 83., 53., 26.],
[ 0., 1., 4., 0.],
[ 0., 0., 0., 0.],
[ 0., 1., 2., 0.]]) / 2, 0.001)
np.testing.assert_allclose(cycle_matrix(data, 4, 4, rainflow_func=rainflow_astm)[0], np.array([[24., 83., 53., 26.],
[0., 1., 4., 0.],
[0., 0., 0., 0.],
[0., 1., 2., 0.]]) / 2, 0.001)
def test_astm_weighted(self):
data = Hawc2io.ReadHawc2(testfilepath + "test").ReadBinary([2]).flatten()
np.testing.assert_allclose(cycle_matrix([(1, data),(1,data)], 4, 4, rainflow_func=rainflow_astm)[0], np.array([[ 24., 83., 53., 26.],
[ 0., 1., 4., 0.],
[ 0., 0., 0., 0.],
[ 0., 1., 2., 0.]]) , 0.001)
np.testing.assert_allclose(cycle_matrix([(1, data), (1, data)], 4, 4, rainflow_func=rainflow_astm)[0], np.array([[24., 83., 53., 26.],
[0., 1.,
4., 0.],
[0., 0.,
0., 0.],
[0., 1., 2., 0.]]), 0.001)
def test_astm_matlab_example(self):
# example from https://se.mathworks.com/help/signal/ref/rainflow.html
fs = 512
X = np.array([-2, 1, -3, 5, -1, 3, -4, 4, -2])
Y = -np.diff(X)[:, np.newaxis] / 2. * np.cos(np.pi *
np.arange(0, 1, 1 / fs))[np.newaxis] + ((X[:-1] + X[1:]) / 2)[:, np.newaxis]
Y = np.r_[Y.flatten(), X[-1]]
range_lst, mean_lst = (rainflow_astm(Y))
np.testing.assert_array_equal(range_lst, [3, 4, 4, 4, 8, 9, 8, 6])
np.testing.assert_array_equal(mean_lst, [-.5, -1, 1, 1, 1, .5, 0, 1])
if 0:
import matplotlib.pyplot as plt
plt.plot(np.arange(0, len(X) - 1 + 1 / fs, 1 / fs), Y)
plt.plot(np.arange(len(X)), X, 'o')
plt.show()
if __name__ == "__main__":
#import sys;sys.argv = ['', 'Test.testName']
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment