added fourier (not finished)

added TurningPoints.rainflow_filter
+translated some parts of chapter4.py
master
Per.Andreas.Brodtkorb 14 years ago
parent 478846ee01
commit 993d888af2

@ -3,7 +3,7 @@
<pydev_project>
<pydev_pathproperty name="org.python.pydev.PROJECT_SOURCE_PATH">
<path>/pywafo/src</path>
<path>/google_pywafo/src</path>
</pydev_pathproperty>
<pydev_property name="org.python.pydev.PYTHON_PROJECT_VERSION">python 2.6</pydev_property>
<pydev_property name="org.python.pydev.PYTHON_PROJECT_INTERPRETER">Default</pydev_property>

@ -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);

@ -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)

@ -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):

Loading…
Cancel
Save