From 993d888af2e6ac740a0d2370e20016020ecac5bf Mon Sep 17 00:00:00 2001 From: "Per.Andreas.Brodtkorb" Date: Mon, 31 Jan 2011 13:30:16 +0000 Subject: [PATCH] added fourier (not finished) added TurningPoints.rainflow_filter +translated some parts of chapter4.py --- pywafo/.pydevproject | 2 +- .../src/wafo/doc/tutorial_scripts/chapter4.py | 199 ++++++++++-------- pywafo/src/wafo/misc.py | 173 ++++++++++++++- pywafo/src/wafo/objects.py | 62 +++++- 4 files changed, 333 insertions(+), 103 deletions(-) diff --git a/pywafo/.pydevproject b/pywafo/.pydevproject index 56f4b02..4884fb8 100644 --- a/pywafo/.pydevproject +++ b/pywafo/.pydevproject @@ -3,7 +3,7 @@ -/pywafo/src +/google_pywafo/src python 2.6 Default diff --git a/pywafo/src/wafo/doc/tutorial_scripts/chapter4.py b/pywafo/src/wafo/doc/tutorial_scripts/chapter4.py index f5f828f..90d84d6 100644 --- a/pywafo/src/wafo/doc/tutorial_scripts/chapter4.py +++ b/pywafo/src/wafo/doc/tutorial_scripts/chapter4.py @@ -1,64 +1,77 @@ -%% CHAPTER4 contains the commands used in Chapter 4 of the tutorial -% -% CALL: Chapter4 -% -% Some of the commands are edited for fast computation. -% Each set of commands is followed by a 'pause' command. -% -% This routine also can print the figures; -% For printing the figures on directory ../bilder/ edit the file and put -% printing=1; - -% Tested on Matlab 5.3 -% History -% Revised pab sept2005 -% Added sections -> easier to evaluate using cellmode evaluation. -% revised pab Feb2004 -% updated call to lc2sdat -% Created by GL July 13, 2000 -% from commands used in Chapter 4 -% +import numpy as np +from scipy import * +from pylab import * + +#! CHAPTER4 contains the commands used in Chapter 4 of the tutorial +#!================================================================= +#! +#! CALL: Chapter4 +#! +#! Some of the commands are edited for fast computation. +#! Each set of commands is followed by a 'pause' command. +#! +#! This routine also can print the figures; +#! For printing the figures on directory ../bilder/ edit the file and put +#! printing=1; + +#! Tested on Matlab 5.3 +#! History +#! Revised pab sept2005 +#! Added sections -> easier to evaluate using cellmode evaluation. +#! revised pab Feb2004 +#! updated call to lc2sdat +#! Created by GL July 13, 2000 +#! from commands used in Chapter 4 +#! -%% Chapter 4 Fatigue load analysis and rain-flow cycles - -pstate = 'off'; +#! Chapter 4 Fatigue load analysis and rain-flow cycles +#!------------------------------------------------------ printing=0; -%set(0,'DefaultAxesFontSize',15') - -%% Section 4.3.1 Crossing intensity -xx_sea = load('sea.dat'); -tp_sea = dat2tp(xx_sea); -lc_sea = tp2lc(tp_sea); -T_sea = xx_sea(end,1)-xx_sea(1,1); -lc_sea(:,2) = lc_sea(:,2)/T_sea; -clf -subplot(221), plot(lc_sea(:,1),lc_sea(:,2)) -title('Crossing intensity, (u, \mu(u))') -subplot(222), semilogx(lc_sea(:,2),lc_sea(:,1)) -title('Crossing intensity, (log \mu(u), u)') -wafostamp([],'(ER)') -disp('Block 1'), pause(pstate) -m_sea = mean(xx_sea(:,2)); -f0_sea = interp1(lc_sea(:,1),lc_sea(:,2),m_sea,'linear') -extr_sea = length(tp_sea)/(2*T_sea); + +#! Section 4.3.1 Crossing intensity +#!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +import wafo.data as wd +import wafo.objects as wo + +xx_sea = wd.sea() +ts = wo.mat2timeseries(xx_sea) +tp = ts.turning_points() +mM = tp.cycle_pairs(kind='min2max') +lc = mM.level_crossings() +lci = lc.copy() +T_sea = lci.args[-1]-lci.args[0] +lci.data = lci.data/T_sea +lci.labels.ylab='Crossing intensity' +subplot(2,2,1) +lci.plot() +subplot(2,2,2) +lci.setplotter(plotmethod='step') +show() + + +m_sea = ts.data.mean() +f0_sea = interp(m_sea, lci.args,lci.data) +extr_sea = len(tp.data)/(2*T_sea) alfa_sea = f0_sea/extr_sea -disp('Block 2'),pause(pstate) +print('alfa = %g ' % alfa_sea ) -%% Section 4.3.2 Extraction of rainflow cycles -%% Min-max and rainflow cycle plots -RFC_sea=tp2rfc(tp_sea); +#! Section 4.3.2 Extraction of rainflow cycles +#!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +#!#! Min-max and rainflow cycle plots + +mM_rfc=tp.cycle_pairs(h=0.3) mM_sea=tp2mm(tp_sea); clf -subplot(122), ccplot(mM_sea); +subplot(122), +mM.plot() title('min-max cycle count') -subplot(121), ccplot(RFC_sea); +subplot(121), +mM_rfc.plot() title('Rainflow cycle count') -wafostamp([],'(ER)') -disp('Block 3'),pause(pstate) -%% Min-max and rainflow cycle distributions +#!#! Min-max and rainflow cycle distributions ampmM_sea=cc2amp(mM_sea); ampRFC_sea=cc2amp(RFC_sea); clf @@ -69,8 +82,8 @@ title('Rainflow amplitude distribution') wafostamp([],'(ER)') disp('Block 4'),pause(pstate) -%% Section 4.3.3 Simulation of rainflow cycles -%% Simulation of cycles in a Markov model +#!#! Section 4.3.3 Simulation of rainflow cycles +#!#! Simulation of cycles in a Markov model n=41; param_m=[-1 1 n]; param_D=[1 n n]; u_markov=levels(param_m); G_markov=mktestmat(param_m,[-0.2 0.2],0.15,1); @@ -84,8 +97,8 @@ wafostamp([],'(ER)') disp('Block 5'),pause(pstate) -%% Rainflow cycles in a transformed Gaussian model -%% Hermite transformed wave data and rainflow filtered turning points, h = 0.2. +#!#! Rainflow cycles in a transformed Gaussian model +#!#! Hermite transformed wave data and rainflow filtered turning points, h = 0.2. me = mean(xx_sea(:,2)); sa = std(xx_sea(:,2)); Hm0_sea = 4*sa; @@ -98,10 +111,10 @@ param_h = [-1.5 2 51]; spec_norm = spec; spec_norm.S = spec_norm.S/sa^2; xx_herm = spec2sdat(spec_norm,[2^15 1],0.1); -% ????? PJ, JR 11-Apr-2001 -% NOTE, in the simulation program spec2sdat -%the spectrum must be normalized to variance 1 -% ????? +#! ????? PJ, JR 11-Apr-2001 +#! NOTE, in the simulation program spec2sdat +#!the spectrum must be normalized to variance 1 +#! ????? h = 0.2; [dtp,u_herm,xx_herm_1]=dat2dtp(param_h,xx_herm,h); clf @@ -112,7 +125,7 @@ title('Rainflow filtered wave data') wafostamp([],'(ER)') disp('Block 6'),pause(pstate) -%% Rainflow cycles and rainflow filtered rainflow cycles in the transformed Gaussian process. +#!#! Rainflow cycles and rainflow filtered rainflow cycles in the transformed Gaussian process. tp_herm=dat2tp(xx_herm); RFC_herm=tp2rfc(tp_herm); mM_herm=tp2mm(tp_herm); @@ -129,7 +142,7 @@ end wafostamp([],'(ER)') disp('Block 7'),pause(pstate) -%% Section 4.3.4 Calculating the rainflow matrix +#!#! Section 4.3.4 Calculating the rainflow matrix Grfc_markov=mctp2rfm({G_markov []}); @@ -139,13 +152,13 @@ subplot(122), cmatplot(u_markov,u_markov,Grfc_markov), axis('square') wafostamp([],'(ER)') disp('Block 8'),pause(pstate) -%% +#!#! clf cmatplot(u_markov,u_markov,{G_markov Grfc_markov},3) wafostamp([],'(ER)') disp('Block 9'),pause(pstate) -%% Min-max-matrix and theoretical rainflow matrix for test Markov sequence. +#!#! Min-max-matrix and theoretical rainflow matrix for test Markov sequence. cmatplot(u_markov,u_markov,{G_markov Grfc_markov},4) subplot(121), axis('square'), title('min2max transition matrix') subplot(122), axis('square'), title('Rainflow matrix') @@ -154,7 +167,7 @@ end wafostamp([],'(ER)') disp('Block 10'),pause(pstate) -%% Observed and theoretical rainflow matrix for test Markov sequence. +#!#! Observed and theoretical rainflow matrix for test Markov sequence. n=length(u_markov); Frfc_markov=dtp2rfm(xxD_markov,n); clf @@ -166,7 +179,7 @@ end wafostamp([],'(ER)') disp('Block 11'),pause(pstate) -%% Smoothed observed and calculated rainflow matrix for test Markov sequence. +#!#! Smoothed observed and calculated rainflow matrix for test Markov sequence. tp_markov=dat2tp(xx_markov); RFC_markov=tp2rfc(tp_markov); h=1; @@ -180,16 +193,16 @@ end wafostamp([],'(ER)') disp('Block 12'),pause(pstate) -%% Rainflow matrix from spectrum +#!#! Rainflow matrix from spectrum clf -%GmM3_herm=spec2mmtpdf(spec,[],'Mm',[],[],2); +#!GmM3_herm=spec2mmtpdf(spec,[],'Mm',[],[],2); GmM3_herm=spec2cmat(spec,[],'Mm',[],param_h,2); pdfplot(GmM3_herm) wafostamp([],'(ER)') disp('Block 13'),pause(pstate) -%% Min-max matrix and theoretical rainflow matrix for Hermite-transformed Gaussian waves. +#!#! Min-max matrix and theoretical rainflow matrix for Hermite-transformed Gaussian waves. Grfc_herm=mctp2rfm({GmM3_herm.f []}); u_herm=levels(param_h); clf @@ -201,7 +214,7 @@ end wafostamp([],'(ER)') disp('Block 14'),pause(pstate) -%% +#!#! clf Grfc_direct_herm=spec2cmat(spec,[],'rfc',[],[],2); subplot(121), pdfplot(GmM3_herm), axis('square'), hold on @@ -212,8 +225,8 @@ wafostamp([],'(ER)') disp('Block 15'),pause(pstate) -%% Observed smoothed and theoretical min-max matrix, -%% (and observed smoothed and theoretical rainflow matrix for Hermite-transformed Gaussian waves). +#!#! Observed smoothed and theoretical min-max matrix, +#!#! (and observed smoothed and theoretical rainflow matrix for Hermite-transformed Gaussian waves). tp_herm=dat2tp(xx_herm); RFC_herm=tp2rfc(tp_herm); mM_herm=tp2mm(tp_herm); @@ -233,10 +246,10 @@ end wafostamp([],'(ER)') disp('Block 16'),pause(pstate) -%% Section 4.3.5 Simulation from crossings and rainflow structure +#!#! Section 4.3.5 Simulation from crossings and rainflow structure -%% Crossing spectrum (smooth curve) and obtained spectrum (wiggled curve) -%% for simulated process with irregularity factor 0.25. +#!#! Crossing spectrum (smooth curve) and obtained spectrum (wiggled curve) +#!#! for simulated process with irregularity factor 0.25. clf cross_herm=dat2lc(xx_herm); alpha1=0.25; @@ -258,8 +271,8 @@ end wafostamp([],'(ER)') disp('Block 16'),pause(pstate) -%% Crossing spectrum (smooth curve) and obtained spectrum (wiggled curve) -%% for simulated process with irregularity factor 0.75. +#!#! Crossing spectrum (smooth curve) and obtained spectrum (wiggled curve) +#!#! for simulated process with irregularity factor 0.75. xx_herm_sim2=lc2sdat(cross_herm,500,alpha2); cross_herm_sim2=dat2lc(xx_herm_sim2); subplot(211) @@ -277,15 +290,15 @@ end wafostamp([],'(ER)') disp('Block 17'),pause(pstate) -%% Section 4.4 Fatigue damage and fatigue life distribution -%% Section 4.4.1 Introduction +#!#! Section 4.4 Fatigue damage and fatigue life distribution +#!#! Section 4.4.1 Introduction beta=3.2; gam=5.5E-10; T_sea=xx_sea(end,1)-xx_sea(1,1); d_beta=cc2dam(RFC_sea,beta)/T_sea; -time_fail=1/gam/d_beta/3600 %in hours of the specific storm +time_fail=1/gam/d_beta/3600 #!in hours of the specific storm disp('Block 18'),pause(pstate) -%% Section 4.4.2 Level crossings -%% Crossing intensity as calculated from the Markov matrix (solid curve) and from the observed rainflow matrix (dashed curve). +#!#! Section 4.4.2 Level crossings +#!#! Crossing intensity as calculated from the Markov matrix (solid curve) and from the observed rainflow matrix (dashed curve). clf mu_markov=cmat2lc(param_m,Grfc_markov); muObs_markov=cmat2lc(param_m,Frfc_markov/(T_markov/2)); @@ -297,8 +310,8 @@ end wafostamp([],'(ER)') disp('Block 19'),pause(pstate) -%% Section 4.4.3 Damage -%% Distribution of damage from different RFC cycles, from calculated theoretical and from observed rainflow matrix. +#!#! Section 4.4.3 Damage +#!#! Distribution of damage from different RFC cycles, from calculated theoretical and from observed rainflow matrix. beta = 4; Dam_markov = cmat2dam(param_m,Grfc_markov,beta) DamObs1_markov = cc2dam(RFC_markov,beta)/(T_markov/2) @@ -318,24 +331,24 @@ wafostamp([],'(ER)') disp('Block 21'),pause(pstate) -%% -%Damplus_markov = lc2dplus(mu_markov,beta) +#!#! +#!Damplus_markov = lc2dplus(mu_markov,beta) pause(pstate) -%% Section 4.4.4 Estimation of S-N curve +#!#! Section 4.4.4 Estimation of S-N curve -%% Load SN-data and plot in log-log scale. +#!#! Load SN-data and plot in log-log scale. SN = load('sn.dat'); s = SN(:,1); N = SN(:,2); clf loglog(N,s,'o'), axis([0 14e5 10 30]) -%if (printing==1), print -deps ../bilder/fatigue_?.eps end +#!if (printing==1), print -deps ../bilder/fatigue_?.eps end wafostamp([],'(ER)') disp('Block 22'),pause(pstate) -%% Check of S-N-model on normal probability paper. +#!#! Check of S-N-model on normal probability paper. normplot(reshape(log(N),8,5)) if (printing==1), print -deps ../bilder/fatigue_17.eps @@ -343,7 +356,7 @@ end wafostamp([],'(ER)') disp('Block 23'),pause(pstate) -%% Estimation of S-N-model on linear scale. +#!#! Estimation of S-N-model on linear scale. clf [e0,beta0,s20] = snplot(s,N,12); title('S-N-data with estimated N(s)','FontSize',20) @@ -353,7 +366,7 @@ end wafostamp([],'(ER)') disp('Block 24'),pause(pstate) -%% Estimation of S-N-model on log-log scale. +#!#! Estimation of S-N-model on log-log scale. clf [e0,beta0,s20] = snplot(s,N,14); title('S-N-data with estimated N(s)','FontSize',20) @@ -363,8 +376,8 @@ end wafostamp([],'(ER)') disp('Block 25'),pause(pstate) -%% Section 4.4.5 From S-N curve to fatigue life distribution -%% Damage intensity as function of $\beta$ +#!#! Section 4.4.5 From S-N curve to fatigue life distribution +#!#! Damage intensity as function of $\beta$ beta = 3:0.1:8; DRFC = cc2dam(RFC_sea,beta); dRFC = DRFC/T_sea; @@ -375,7 +388,7 @@ end wafostamp([],'(ER)') disp('Block 26'),pause(pstate) -%% Fatigue life distribution with sea load. +#!#! Fatigue life distribution with sea load. dam0 = cc2dam(RFC_sea,beta0)/T_sea; [t0,F0] = ftf(e0,dam0,s20,0.5,1); [t1,F1] = ftf(e0,dam0,s20,0,1); diff --git a/pywafo/src/wafo/misc.py b/pywafo/src/wafo/misc.py index 8092d75..62751a8 100644 --- a/pywafo/src/wafo/misc.py +++ b/pywafo/src/wafo/misc.py @@ -1,12 +1,12 @@ ''' -Misc lsdkfalsdflasdfl +Misc ''' from __future__ import division import sys import fractions import numpy as np -from numpy import (abs, amax, any, logical_and, arange, linspace, atleast_1d, +from numpy import (abs, amax, any, logical_and, arange, linspace, atleast_1d,atleast_2d, array, asarray, broadcast_arrays, ceil, floor, frexp, hypot, sqrt, arctan2, sin, cos, exp, log, mod, diff, empty_like, finfo, inf, pi, interp, isnan, isscalar, zeros, ones, @@ -16,6 +16,7 @@ import types import warnings from wafo import plotbackend + try: import wafo.c_library as clib except: @@ -31,6 +32,7 @@ __all__ = ['is_numlike','JITImport', 'DotDict', 'Bunch', 'printf', 'sub_dict_sel 'discretize', 'discretize2', 'pol2cart', 'cart2pol', 'meshgrid', 'ndgrid', 'trangood', 'tranproc', 'histgrm', 'num2pistr'] + def is_numlike(obj): 'return true if *obj* looks like a number' try: @@ -2039,6 +2041,173 @@ def num2pistr(x, n=3): xtxt = format % x return xtxt +def fourier(x,t,T=None,M=None,N=None, method='trapz'): + ''' + Returns Fourier coefficients. + + CALL: [a,b] = fourier(t,x,T,M); + + + Parameters + ---------- + x : array-like + vector or matrix of row vectors with data points shape p x n. + t : array-like + vector with n values indexed from 1 to N. + T : real scalar + primitive period of signal, i.e., smallest period. (default T = t[-1]-t[0] + M : scalar integer + defines no of harmonics desired (default M = N) + N : scalar integer + no of data points (default len(t)) + + Returns + ------- + a,b = Fourier coefficients size M x P + + FOURIER finds the coefficients for a Fourier series representation + of the signal x(t) (given in digital form). It is assumed the signal + is periodic over T. N is the number of data points, and M-1 is the + number of coefficients. + + The signal can be estimated by using M-1 harmonics by: + M-1 + x[i] = 0.5*a[0] + sum (a[n]*c[n,i] + b[n]*s[n,i]) + n=1 + where + c[n,i] = cos(2*pi*(n-1)*t[i]/T) + s[n,i] = sin(2*pi*(n-1)*t[i]/T) + + Note that a[0] is the "dc value". + Remaining values are a[1], a[2], ... , a[M-1]. + + Example + ------- + >>> t = linspace(0,4*T) + >>> x = sin(t); + >>> a, b = fourier(t, x, T=2*pi, M=5) + + See also + -------- + fft + ''' + x = np.atleast_2d(x) + p, n = x.shape + if t is None: + t = np.arange(n) + else: + t = np.atleast_1d(t) + + n = len(t) if n is None else n + m = n if n is none else m + T = t[-1]-t[0] if T is None else T + + + +# switch 0 +# case -1, +# % Define the vectors for computing the Fourier coefficients +# % +# a = zeros(M,P); +# b = zeros(M,P); +# a(1,:) = simpson(x); +# +# % +# % Compute M-1 more coefficients +# tmp = 2*pi*t(:,ones(1,P))/T; +# % tmp = 2*pi*(0:N-1).'/(N-1); +# for n1 = 1:M-1, +# n = n1+1; +# a(n,:) = simpson(x.*cos(n1*tmp)); +# b(n,:) = simpson(x.*sin(n1*tmp)); +# end +# +# a = 2*a/N; +# b = 2*b/N; +# +# case 0, +# % +# a = zeros(M,P); +# b = zeros(M,P); +# a(1,:) = trapz(t,x); +# +# % +# % Compute M-1 more coefficients +# tmp = 2*pi*t(:,ones(1,P))/T; +# % tmp = 2*pi*(0:N-1).'/(N-1); +# for n1 = 1:M-1, +# n = n1+1; +# a(n,:) = trapz(t,x.*cos(n1*tmp)); +# b(n,:) = trapz(t,x.*sin(n1*tmp)); +# end +# a = a/pi; +# b = b/pi; +# +# case 1, +# % Define the vectors for computing the Fourier coefficients +# % +# a = zeros(M,P); +# b = zeros(M,P); +# +# % +# % Compute the dc-level (the a(0) component). +# % +# % Note: the index has to begin with "1". +# % +# +# a(1,:) = sum(x); +# +# % +# % Compute M-1 more coefficients +# tmp = 2*pi*t(:,ones(1,P))/T; +# % tmp = 2*pi*(0:N-1).'/(N-1); +# for n1 = 1:M-1, +# n = n1+1; +# a(n,:) = sum(x.*cos(n1*tmp)); +# b(n,:) = sum(x.*sin(n1*tmp)); +# end +# a = 2*a/N; +# b = 2*b/N; +# case 2, +# % Define the vectors for computing the Fourier coefficients +# % +# a = zeros(M,P); +# b = zeros(M,P); +# a(1,:) = trapz(x); +# +# % +# % Compute M-1 more coefficients +# tmp = 2*pi*t(:,ones(1,P))/T; +# % tmp = 2*pi*(0:N-1).'/(N-1); +# for n1 = 1:M-1, +# n = n1+1; +# a(n,:) = trapz(x.*cos(n1*tmp)); +# b(n,:) = trapz(x.*sin(n1*tmp)); +# end +# +# a = 2*a/N; +# b = 2*b/N; +# case 3 +# % Alternative: faster for large M, but gives different results than above. +# nper = diff(t([1 end]))/T; %No of periods given +# if nper == round(nper), +# N1 = N/nper; +# else +# N1 = N; +# end +# +# % Fourier coefficients by fft +# Fcof1 = 2*ifft(x(1:N1,:),[],1); +# Pcor = [1; exp(sqrt(-1)*(1:M-1).'*t(1))]; % correction term to get +# % the correct integration limits +# Fcof = Fcof1(1:M,:).*Pcor(:,ones(1,P)); +# a = real(Fcof(1:M,:)); +# b = imag(Fcof(1:M,:)); +# end +# return + + + def _test_find_cross(): t = findcross([0, 0, 1, -1, 1], 0) diff --git a/pywafo/src/wafo/objects.py b/pywafo/src/wafo/objects.py index 14c1f50..f46aade 100644 --- a/pywafo/src/wafo/objects.py +++ b/pywafo/src/wafo/objects.py @@ -16,7 +16,7 @@ from __future__ import division from wafo.transform.core import TrData from wafo.transform.models import TrHermite, TrOchi, TrLinear from wafo.stats import edf -from wafo.misc import (nextpow2, findtp, findtc, findcross, +from wafo.misc import (nextpow2, findtp, findrfc, findtc, findcross, ecross, JITImport, DotDict, gravity) from wafodata import WafoData from wafo.interpolate import SmoothSpline @@ -632,7 +632,50 @@ class TurningPoints(WafoData): self.args = ravel(self.args) self.data = ravel(self.data) - def cycle_pairs(self, kind='min2max'): + def rainflow_filter(self, h=0.0, method='clib'): + ''' + Return rainflow filtered turning points (tp). + + Parameters + ---------- + h : scalar + a threshold + if h<=0, then tp is a sequence of turning points (default) + if h>0, then all rainflow cycles with height smaller than + h are removed. + + Returns + ------- + tp : TurningPoints object + with times and turning points. + + Example: + >>> import wafo.data + >>> x = wafo.data.sea() + >>> x1 = x[:200,:] + >>> ts1 = mat2timeseries(x1) + >>> tp = ts1.turning_points(wavetype='Mw') + >>> tph = tp.rainflow_filter(h=0.3) + >>> hs = ts1.plot() + >>> hp = tp.plot('ro') + >>> hph = tph.plot('k.') + + See also + --------- + findcross, + findrfc + findtp + ''' + ind = findrfc(self.data, max(h, 0.0), method) + try: + t = self.args[ind] + except: + t = ind + mean = self.mean() + sigma = self.sigma + return TurningPoints(self.data[ind], t, mean=mean, sigma=sigma) + + def cycle_pairs(self, h=0, kind='min2max', method='clib'): """ Return min2Max or Max2min cycle pairs from turning points Parameters @@ -660,7 +703,12 @@ class TurningPoints(WafoData): TurningPoints SurvivalCycleCount """ - if self.data[0] > self.data[1]: + if h>0: + ind = findrfc(self.data, h, method=method) + data = self.data(ind) + else: + data = self.data + if data[0] > data[1]: im = 1 iM = 0 else: @@ -670,12 +718,12 @@ class TurningPoints(WafoData): # Extract min-max and max-min cycle pairs #n = len(self.data) if kind.lower().startswith('min2max'): - m = self.data[im:-1:2] - M = self.data[im + 1::2] + m = data[im:-1:2] + M = data[im + 1::2] else: kind = 'max2min' - M = self.data[iM:-1:2] - m = self.data[iM + 1::2] + M = data[iM:-1:2] + m = data[iM + 1::2] return CyclePairs(M, m, kind=kind, mean=self.mean, sigma=self.sigma) def mat2timeseries(x):