danh
/
pywafo
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

misc fixes

Browse Source
master
Per.Andreas.Brodtkorb 14 years ago
parent 993d888af2
commit 62316ed499
5 changed files with 508 additions and 497 deletions
Show Stats Download Patch File Download Diff File

3
pywafo/src/wafo/doc/tutorial_scripts/chapter1.py
Unescape Escape View File

@ -64,7 +64,8 @@ dtyex = S1.to_t_pdf(pdef='Tt', paramt=(0, 10, 51), nit=3)
dtyest = Sest.to_t_pdf(pdef='Tt', paramt=(0, 10, 51), nit=3) dtyest = Sest.to_t_pdf(pdef='Tt', paramt=(0, 10, 51), nit=3)
T, index = ts.wave_periods(vh=0, pdef='d2u') T, index = ts.wave_periods(vh=0, pdef='d2u')
wm.histgrm(T, n=25, odd=True, scale=True) bins = wm.good_bins(T, num_bins=25, odd=True)
wm.plot_histgrm(T, bins=bins, scale=True)
dtyex.plot() dtyex.plot()
dtyest.plot('-.') dtyest.plot('-.')

2
pywafo/src/wafo/doc/tutorial_scripts/chapter3.py
Unescape Escape View File

@ -37,7 +37,7 @@ t = linspace(0.01,8,200);
ftc = wk.TKDE(Tc, L2=0, inc=128) ftc = wk.TKDE(Tc, L2=0, inc=128)
plot(t,ftc.eval_grid(t), t, ftc.eval_grid_fast(t),'-.') plot(t,ftc.eval_grid(t), t, ftc.eval_grid_fast(t),'-.')
wm.histgrm(Tc,scale=True) wm.plot_histgrm(Tc,normed=True)
title('Kernel Density Estimates') title('Kernel Density Estimates')
axis([0, 8, 0, 0.5]) axis([0, 8, 0, 0.5])
show() show()

675
pywafo/src/wafo/doc/tutorial_scripts/chapter4.py
Unescape Escape View File

@ -39,363 +39,368 @@ xx_sea = wd.sea()
ts = wo.mat2timeseries(xx_sea) ts = wo.mat2timeseries(xx_sea)
tp = ts.turning_points() tp = ts.turning_points()
mM = tp.cycle_pairs(kind='min2max') mM = tp.cycle_pairs(kind='min2max')
lc = mM.level_crossings() lc = mM.level_crossings(intensity=True)
lci = lc.copy() T_sea = ts.args[-1]-ts.args[0]
T_sea = lci.args[-1]-lci.args[0]
lci.data = lci.data/T_sea subplot(1,2,1)
lci.labels.ylab='Crossing intensity' lc.plot()
subplot(2,2,1) subplot(1,2,2)
lci.plot() lc.setplotter(plotmethod='step')
subplot(2,2,2) lc.plot()
lci.setplotter(plotmethod='step')
show() show()
m_sea = ts.data.mean() m_sea = ts.data.mean()
f0_sea = interp(m_sea, lci.args,lci.data) f0_sea = interp(m_sea, lc.args,lc.data)
extr_sea = len(tp.data)/(2*T_sea) extr_sea = len(tp.data)/(2*T_sea)
alfa_sea = f0_sea/extr_sea alfa_sea = f0_sea/extr_sea
print('alfa = %g ' % alfa_sea ) print('alfa = %g ' % alfa_sea )
#! Section 4.3.2 Extraction of rainflow cycles #! Section 4.3.2 Extraction of rainflow cycles
#!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#!#! Min-max and rainflow cycle plots #! Min-max and rainflow cycle plots
#!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
mM_rfc = tp.cycle_pairs(h=0.3)
mM_rfc=tp.cycle_pairs(h=0.3) clf()
mM_sea=tp2mm(tp_sea);
clf
subplot(122), subplot(122),
mM.plot() mM.plot()
title('min-max cycle count') title('min-max cycle count')
subplot(121), subplot(121),
mM_rfc.plot() mM_rfc.plot()
title('Rainflow cycle count') title('Rainflow cycle count')
show()
#!#! Min-max and rainflow cycle distributions #! Min-max and rainflow cycle distributions
ampmM_sea=cc2amp(mM_sea); #!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
ampRFC_sea=cc2amp(RFC_sea); import wafo.misc as wm
clf ampmM_sea = mM.amplitudes()
subplot(221), hist(ampmM_sea,25); ampRFC_sea = mM_rfc.amplitudes()
clf()
subplot(121)
wm.plot_histgrm(ampmM_sea,25)
ylim = gca().get_ylim()
title('min-max amplitude distribution') title('min-max amplitude distribution')
subplot(222), hist(ampRFC_sea,25); subplot(122)
wm.plot_histgrm(ampRFC_sea,25)
gca().set_ylim(ylim)
title('Rainflow amplitude distribution') title('Rainflow amplitude distribution')
wafostamp([],'(ER)') show()
disp('Block 4'),pause(pstate)
#!#! Section 4.3.3 Simulation of rainflow cycles #!#! Section 4.3.3 Simulation of rainflow cycles
#!#! Simulation of cycles in a Markov model #!#! Simulation of cycles in a Markov model
n=41; param_m=[-1 1 n]; param_D=[1 n n]; n=41; param_m=[-1, 1, n]; param_D=[1, n, n];
u_markov=levels(param_m); u_markov=levels(param_m);
G_markov=mktestmat(param_m,[-0.2 0.2],0.15,1); G_markov=mktestmat(param_m,[-0.2, 0.2],0.15,1);
T_markov=5000; T_markov=5000;
xxD_markov=mctpsim({G_markov []},T_markov); #xxD_markov=mctpsim({G_markov [,]},T_markov);
xx_markov=[(1:T_markov)' u_markov(xxD_markov)']; #xx_markov=[(1:T_markov)' u_markov(xxD_markov)'];
clf #clf
plot(xx_markov(1:50,1),xx_markov(1:50,2)) #plot(xx_markov(1:50,1),xx_markov(1:50,2))
title('Markov chain of turning points') #title('Markov chain of turning points')
wafostamp([],'(ER)') #wafostamp([],'(ER)')
disp('Block 5'),pause(pstate) #disp('Block 5'),pause(pstate)
#
#
#!#! Rainflow cycles in a transformed Gaussian model ##!#! Rainflow cycles in a transformed Gaussian model
#!#! Hermite transformed wave data and rainflow filtered turning points, h = 0.2. ##!#! Hermite transformed wave data and rainflow filtered turning points, h = 0.2.
me = mean(xx_sea(:,2)); #me = mean(xx_sea(:,2));
sa = std(xx_sea(:,2)); #sa = std(xx_sea(:,2));
Hm0_sea = 4*sa; #Hm0_sea = 4*sa;
Tp_sea = 1/max(lc_sea(:,2)); #Tp_sea = 1/max(lc_sea(:,2));
spec = jonswap([],[Hm0_sea Tp_sea]); #spec = jonswap([],[Hm0_sea Tp_sea]);
#
[sk, ku] = spec2skew(spec); #[sk, ku] = spec2skew(spec);
spec.tr = hermitetr([],[sa sk ku me]); #spec.tr = hermitetr([],[sa sk ku me]);
param_h = [-1.5 2 51]; #param_h = [-1.5 2 51];
spec_norm = spec; #spec_norm = spec;
spec_norm.S = spec_norm.S/sa^2; #spec_norm.S = spec_norm.S/sa^2;
xx_herm = spec2sdat(spec_norm,[2^15 1],0.1); #xx_herm = spec2sdat(spec_norm,[2^15 1],0.1);
#! ????? PJ, JR 11-Apr-2001 ##! ????? PJ, JR 11-Apr-2001
#! NOTE, in the simulation program spec2sdat ##! NOTE, in the simulation program spec2sdat
#!the spectrum must be normalized to variance 1 ##!the spectrum must be normalized to variance 1
#! ????? ##! ?????
h = 0.2; #h = 0.2;
[dtp,u_herm,xx_herm_1]=dat2dtp(param_h,xx_herm,h); #[dtp,u_herm,xx_herm_1]=dat2dtp(param_h,xx_herm,h);
clf #clf
plot(xx_herm(:,1),xx_herm(:,2),'k','LineWidth',2); hold on; #plot(xx_herm(:,1),xx_herm(:,2),'k','LineWidth',2); hold on;
plot(xx_herm_1(:,1),xx_herm_1(:,2),'k--','Linewidth',2); #plot(xx_herm_1(:,1),xx_herm_1(:,2),'k--','Linewidth',2);
axis([0 50 -1 1]), hold off; #axis([0 50 -1 1]), hold off;
title('Rainflow filtered wave data') #title('Rainflow filtered wave data')
wafostamp([],'(ER)') #wafostamp([],'(ER)')
disp('Block 6'),pause(pstate) #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); #tp_herm=dat2tp(xx_herm);
RFC_herm=tp2rfc(tp_herm); #RFC_herm=tp2rfc(tp_herm);
mM_herm=tp2mm(tp_herm); #mM_herm=tp2mm(tp_herm);
h=0.2; #h=0.2;
[dtp,u,tp_herm_1]=dat2dtp(param_h,xx_herm,h); #[dtp,u,tp_herm_1]=dat2dtp(param_h,xx_herm,h);
RFC_herm_1 = tp2rfc(tp_herm_1); #RFC_herm_1 = tp2rfc(tp_herm_1);
clf #clf
subplot(121), ccplot(RFC_herm) #subplot(121), ccplot(RFC_herm)
title('h=0') #title('h=0')
subplot(122), ccplot(RFC_herm_1) #subplot(122), ccplot(RFC_herm_1)
title('h=0.2') #title('h=0.2')
if (printing==1), print -deps ../bilder/fatigue_8.eps #if (printing==1), print -deps ../bilder/fatigue_8.eps
end #end
wafostamp([],'(ER)') #wafostamp([],'(ER)')
disp('Block 7'),pause(pstate) #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 []}); #Grfc_markov=mctp2rfm({G_markov []});
clf #clf
subplot(121), cmatplot(u_markov,u_markov,G_markov), axis('square') #subplot(121), cmatplot(u_markov,u_markov,G_markov), axis('square')
subplot(122), cmatplot(u_markov,u_markov,Grfc_markov), axis('square') #subplot(122), cmatplot(u_markov,u_markov,Grfc_markov), axis('square')
wafostamp([],'(ER)') #wafostamp([],'(ER)')
disp('Block 8'),pause(pstate) #disp('Block 8'),pause(pstate)
#
#!#! ##!#!
clf #clf
cmatplot(u_markov,u_markov,{G_markov Grfc_markov},3) #cmatplot(u_markov,u_markov,{G_markov Grfc_markov},3)
wafostamp([],'(ER)') #wafostamp([],'(ER)')
disp('Block 9'),pause(pstate) #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) #cmatplot(u_markov,u_markov,{G_markov Grfc_markov},4)
subplot(121), axis('square'), title('min2max transition matrix') #subplot(121), axis('square'), title('min2max transition matrix')
subplot(122), axis('square'), title('Rainflow matrix') #subplot(122), axis('square'), title('Rainflow matrix')
if (printing==1), print -deps ../bilder/fatigue_9.eps #if (printing==1), print -deps ../bilder/fatigue_9.eps
end #end
wafostamp([],'(ER)') #wafostamp([],'(ER)')
disp('Block 10'),pause(pstate) #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); #n=length(u_markov);
Frfc_markov=dtp2rfm(xxD_markov,n); #Frfc_markov=dtp2rfm(xxD_markov,n);
clf #clf
cmatplot(u_markov,u_markov,{Frfc_markov Grfc_markov*T_markov/2},3) #cmatplot(u_markov,u_markov,{Frfc_markov Grfc_markov*T_markov/2},3)
subplot(121), axis('square'), title('Observed rainflow matrix') #subplot(121), axis('square'), title('Observed rainflow matrix')
subplot(122), axis('square'), title('Theoretical rainflow matrix') #subplot(122), axis('square'), title('Theoretical rainflow matrix')
if (printing==1), print -deps ../bilder/fatigue_10.eps #if (printing==1), print -deps ../bilder/fatigue_10.eps
end #end
wafostamp([],'(ER)') #wafostamp([],'(ER)')
disp('Block 11'),pause(pstate) #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); #tp_markov=dat2tp(xx_markov);
RFC_markov=tp2rfc(tp_markov); #RFC_markov=tp2rfc(tp_markov);
h=1; #h=1;
Frfc_markov_smooth=cc2cmat(param_m,RFC_markov,[],1,h); #Frfc_markov_smooth=cc2cmat(param_m,RFC_markov,[],1,h);
clf #clf
cmatplot(u_markov,u_markov,{Frfc_markov_smooth Grfc_markov*T_markov/2},4) #cmatplot(u_markov,u_markov,{Frfc_markov_smooth Grfc_markov*T_markov/2},4)
subplot(121), axis('square'), title('Smoothed observed rainflow matrix') #subplot(121), axis('square'), title('Smoothed observed rainflow matrix')
subplot(122), axis('square'), title('Theoretical rainflow matrix') #subplot(122), axis('square'), title('Theoretical rainflow matrix')
if (printing==1), print -deps ../bilder/fatigue_11.eps #if (printing==1), print -deps ../bilder/fatigue_11.eps
end #end
wafostamp([],'(ER)') #wafostamp([],'(ER)')
disp('Block 12'),pause(pstate) #disp('Block 12'),pause(pstate)
#
#!#! Rainflow matrix from spectrum ##!#! Rainflow matrix from spectrum
clf #clf
#!GmM3_herm=spec2mmtpdf(spec,[],'Mm',[],[],2); ##!GmM3_herm=spec2mmtpdf(spec,[],'Mm',[],[],2);
GmM3_herm=spec2cmat(spec,[],'Mm',[],param_h,2); #GmM3_herm=spec2cmat(spec,[],'Mm',[],param_h,2);
pdfplot(GmM3_herm) #pdfplot(GmM3_herm)
wafostamp([],'(ER)') #wafostamp([],'(ER)')
disp('Block 13'),pause(pstate) #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 []}); #Grfc_herm=mctp2rfm({GmM3_herm.f []});
u_herm=levels(param_h); #u_herm=levels(param_h);
clf #clf
cmatplot(u_herm,u_herm,{GmM3_herm.f Grfc_herm},4) #cmatplot(u_herm,u_herm,{GmM3_herm.f Grfc_herm},4)
subplot(121), axis('square'), title('min-max matrix') #subplot(121), axis('square'), title('min-max matrix')
subplot(122), axis('square'), title('Theoretical rainflow matrix') #subplot(122), axis('square'), title('Theoretical rainflow matrix')
if (printing==1), print -deps ../bilder/fatigue_12.eps #if (printing==1), print -deps ../bilder/fatigue_12.eps
end #end
wafostamp([],'(ER)') #wafostamp([],'(ER)')
disp('Block 14'),pause(pstate) #disp('Block 14'),pause(pstate)
#
#!#! ##!#!
clf #clf
Grfc_direct_herm=spec2cmat(spec,[],'rfc',[],[],2); #Grfc_direct_herm=spec2cmat(spec,[],'rfc',[],[],2);
subplot(121), pdfplot(GmM3_herm), axis('square'), hold on #subplot(121), pdfplot(GmM3_herm), axis('square'), hold on
subplot(122), pdfplot(Grfc_direct_herm), axis('square'), hold off #subplot(122), pdfplot(Grfc_direct_herm), axis('square'), hold off
if (printing==1), print -deps ../bilder/fig_mmrfcjfr.eps #if (printing==1), print -deps ../bilder/fig_mmrfcjfr.eps
end #end
wafostamp([],'(ER)') #wafostamp([],'(ER)')
disp('Block 15'),pause(pstate) #disp('Block 15'),pause(pstate)
#
#
#!#! Observed smoothed and theoretical min-max matrix, ##!#! Observed smoothed and theoretical min-max matrix,
#!#! (and observed smoothed and theoretical rainflow matrix for Hermite-transformed Gaussian waves). ##!#! (and observed smoothed and theoretical rainflow matrix for Hermite-transformed Gaussian waves).
tp_herm=dat2tp(xx_herm); #tp_herm=dat2tp(xx_herm);
RFC_herm=tp2rfc(tp_herm); #RFC_herm=tp2rfc(tp_herm);
mM_herm=tp2mm(tp_herm); #mM_herm=tp2mm(tp_herm);
h=0.2; #h=0.2;
FmM_herm_smooth=cc2cmat(param_h,mM_herm,[],1,h); #FmM_herm_smooth=cc2cmat(param_h,mM_herm,[],1,h);
Frfc_herm_smooth=cc2cmat(param_h,RFC_herm,[],1,h); #Frfc_herm_smooth=cc2cmat(param_h,RFC_herm,[],1,h);
T_herm=xx_herm(end,1)-xx_herm(1,1); #T_herm=xx_herm(end,1)-xx_herm(1,1);
clf #clf
cmatplot(u_herm,u_herm,{FmM_herm_smooth GmM3_herm.f*length(mM_herm) ; ... #cmatplot(u_herm,u_herm,{FmM_herm_smooth GmM3_herm.f*length(mM_herm) ; ...
Frfc_herm_smooth Grfc_herm*length(RFC_herm)},4) # Frfc_herm_smooth Grfc_herm*length(RFC_herm)},4)
subplot(221), axis('square'), title('Observed smoothed min-max matrix') #subplot(221), axis('square'), title('Observed smoothed min-max matrix')
subplot(222), axis('square'), title('Theoretical min-max matrix') #subplot(222), axis('square'), title('Theoretical min-max matrix')
subplot(223), axis('square'), title('Observed smoothed rainflow matrix') #subplot(223), axis('square'), title('Observed smoothed rainflow matrix')
subplot(224), axis('square'), title('Theoretical rainflow matrix') #subplot(224), axis('square'), title('Theoretical rainflow matrix')
if (printing==1), print -deps ../bilder/fatigue_13.eps #if (printing==1), print -deps ../bilder/fatigue_13.eps
end #end
wafostamp([],'(ER)') #wafostamp([],'(ER)')
disp('Block 16'),pause(pstate) #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) ##!#! Crossing spectrum (smooth curve) and obtained spectrum (wiggled curve)
#!#! for simulated process with irregularity factor 0.25. ##!#! for simulated process with irregularity factor 0.25.
clf #clf
cross_herm=dat2lc(xx_herm); #cross_herm=dat2lc(xx_herm);
alpha1=0.25; #alpha1=0.25;
alpha2=0.75; #alpha2=0.75;
xx_herm_sim1=lc2sdat(cross_herm,500,alpha1); #xx_herm_sim1=lc2sdat(cross_herm,500,alpha1);
cross_herm_sim1=dat2lc(xx_herm_sim1); #cross_herm_sim1=dat2lc(xx_herm_sim1);
subplot(211) #subplot(211)
plot(cross_herm(:,1),cross_herm(:,2)/max(cross_herm(:,2))) #plot(cross_herm(:,1),cross_herm(:,2)/max(cross_herm(:,2)))
hold on #hold on
stairs(cross_herm_sim1(:,1),... #stairs(cross_herm_sim1(:,1),...
cross_herm_sim1(:,2)/max(cross_herm_sim1(:,2))) # cross_herm_sim1(:,2)/max(cross_herm_sim1(:,2)))
hold off #hold off
title('Crossing intensity, \alpha = 0.25') #title('Crossing intensity, \alpha = 0.25')
subplot(212) #subplot(212)
plot(xx_herm_sim1(:,1),xx_herm_sim1(:,2)) #plot(xx_herm_sim1(:,1),xx_herm_sim1(:,2))
title('Simulated load, \alpha = 0.25') #title('Simulated load, \alpha = 0.25')
if (printing==1), print -deps ../bilder/fatigue_14_25.eps #if (printing==1), print -deps ../bilder/fatigue_14_25.eps
end #end
wafostamp([],'(ER)') #wafostamp([],'(ER)')
disp('Block 16'),pause(pstate) #disp('Block 16'),pause(pstate)
#
#!#! Crossing spectrum (smooth curve) and obtained spectrum (wiggled curve) ##!#! Crossing spectrum (smooth curve) and obtained spectrum (wiggled curve)
#!#! for simulated process with irregularity factor 0.75. ##!#! for simulated process with irregularity factor 0.75.
xx_herm_sim2=lc2sdat(cross_herm,500,alpha2); #xx_herm_sim2=lc2sdat(cross_herm,500,alpha2);
cross_herm_sim2=dat2lc(xx_herm_sim2); #cross_herm_sim2=dat2lc(xx_herm_sim2);
subplot(211) #subplot(211)
plot(cross_herm(:,1),cross_herm(:,2)/max(cross_herm(:,2))) #plot(cross_herm(:,1),cross_herm(:,2)/max(cross_herm(:,2)))
hold on #hold on
stairs(cross_herm_sim2(:,1),... #stairs(cross_herm_sim2(:,1),...
cross_herm_sim2(:,2)/max(cross_herm_sim2(:,2))) # cross_herm_sim2(:,2)/max(cross_herm_sim2(:,2)))
hold off #hold off
title('Crossing intensity, \alpha = 0.75') #title('Crossing intensity, \alpha = 0.75')
subplot(212) #subplot(212)
plot(xx_herm_sim2(:,1),xx_herm_sim2(:,2)) #plot(xx_herm_sim2(:,1),xx_herm_sim2(:,2))
title('Simulated load, \alpha = 0.75') #title('Simulated load, \alpha = 0.75')
if (printing==1), print -deps ../bilder/fatigue_14_75.eps #if (printing==1), print -deps ../bilder/fatigue_14_75.eps
end #end
wafostamp([],'(ER)') #wafostamp([],'(ER)')
disp('Block 17'),pause(pstate) #disp('Block 17'),pause(pstate)
#
#!#! Section 4.4 Fatigue damage and fatigue life distribution ##!#! Section 4.4 Fatigue damage and fatigue life distribution
#!#! Section 4.4.1 Introduction ##!#! Section 4.4.1 Introduction
beta=3.2; gam=5.5E-10; T_sea=xx_sea(end,1)-xx_sea(1,1); #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; #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) #disp('Block 18'),pause(pstate)
#
#!#! Section 4.4.2 Level crossings ##!#! Section 4.4.2 Level crossings
#!#! Crossing intensity as calculated from the Markov matrix (solid curve) and from the observed rainflow matrix (dashed curve). ##!#! Crossing intensity as calculated from the Markov matrix (solid curve) and from the observed rainflow matrix (dashed curve).
clf #clf
mu_markov=cmat2lc(param_m,Grfc_markov); #mu_markov=cmat2lc(param_m,Grfc_markov);
muObs_markov=cmat2lc(param_m,Frfc_markov/(T_markov/2)); #muObs_markov=cmat2lc(param_m,Frfc_markov/(T_markov/2));
clf #clf
plot(mu_markov(:,1),mu_markov(:,2),muObs_markov(:,1),muObs_markov(:,2),'--') #plot(mu_markov(:,1),mu_markov(:,2),muObs_markov(:,1),muObs_markov(:,2),'--')
title('Theoretical and observed crossing intensity ') #title('Theoretical and observed crossing intensity ')
if (printing==1), print -deps ../bilder/fatigue_15.eps #if (printing==1), print -deps ../bilder/fatigue_15.eps
end #end
wafostamp([],'(ER)') #wafostamp([],'(ER)')
disp('Block 19'),pause(pstate) #disp('Block 19'),pause(pstate)
#
#!#! Section 4.4.3 Damage ##!#! Section 4.4.3 Damage
#!#! Distribution of damage from different RFC cycles, from calculated theoretical and from observed rainflow matrix. ##!#! Distribution of damage from different RFC cycles, from calculated theoretical and from observed rainflow matrix.
beta = 4; #beta = 4;
Dam_markov = cmat2dam(param_m,Grfc_markov,beta) #Dam_markov = cmat2dam(param_m,Grfc_markov,beta)
DamObs1_markov = cc2dam(RFC_markov,beta)/(T_markov/2) #DamObs1_markov = cc2dam(RFC_markov,beta)/(T_markov/2)
DamObs2_markov = cmat2dam(param_m,Frfc_markov,beta)/(T_markov/2) #DamObs2_markov = cmat2dam(param_m,Frfc_markov,beta)/(T_markov/2)
disp('Block 20'),pause(pstate) #disp('Block 20'),pause(pstate)
#
Dmat_markov = cmat2dmat(param_m,Grfc_markov,beta); #Dmat_markov = cmat2dmat(param_m,Grfc_markov,beta);
DmatObs_markov = cmat2dmat(param_m,Frfc_markov,beta)/(T_markov/2); #DmatObs_markov = cmat2dmat(param_m,Frfc_markov,beta)/(T_markov/2);
clf #clf
subplot(121), cmatplot(u_markov,u_markov,Dmat_markov,4) #subplot(121), cmatplot(u_markov,u_markov,Dmat_markov,4)
title('Theoretical damage matrix') #title('Theoretical damage matrix')
subplot(122), cmatplot(u_markov,u_markov,DmatObs_markov,4) #subplot(122), cmatplot(u_markov,u_markov,DmatObs_markov,4)
title('Observed damage matrix') #title('Observed damage matrix')
if (printing==1), print -deps ../bilder/fatigue_16.eps #if (printing==1), print -deps ../bilder/fatigue_16.eps
end #end
wafostamp([],'(ER)') #wafostamp([],'(ER)')
disp('Block 21'),pause(pstate) #disp('Block 21'),pause(pstate)
#
#
#!#! ##!#!
#!Damplus_markov = lc2dplus(mu_markov,beta) ##!Damplus_markov = lc2dplus(mu_markov,beta)
pause(pstate) #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'); #SN = load('sn.dat');
s = SN(:,1); #s = SN(:,1);
N = SN(:,2); #N = SN(:,2);
clf #clf
loglog(N,s,'o'), axis([0 14e5 10 30]) #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)') #wafostamp([],'(ER)')
disp('Block 22'),pause(pstate) #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)) #normplot(reshape(log(N),8,5))
if (printing==1), print -deps ../bilder/fatigue_17.eps #if (printing==1), print -deps ../bilder/fatigue_17.eps
end #end
wafostamp([],'(ER)') #wafostamp([],'(ER)')
disp('Block 23'),pause(pstate) #disp('Block 23'),pause(pstate)
#
#!#! Estimation of S-N-model on linear scale. ##!#! Estimation of S-N-model on linear scale.
clf #clf
[e0,beta0,s20] = snplot(s,N,12); #[e0,beta0,s20] = snplot(s,N,12);
title('S-N-data with estimated N(s)','FontSize',20) #title('S-N-data with estimated N(s)','FontSize',20)
set(gca,'FontSize',20) #set(gca,'FontSize',20)
if (printing==1), print -deps ../bilder/fatigue_18a.eps #if (printing==1), print -deps ../bilder/fatigue_18a.eps
end #end
wafostamp([],'(ER)') #wafostamp([],'(ER)')
disp('Block 24'),pause(pstate) #disp('Block 24'),pause(pstate)
#
#!#! Estimation of S-N-model on log-log scale. ##!#! Estimation of S-N-model on log-log scale.
clf #clf
[e0,beta0,s20] = snplot(s,N,14); #[e0,beta0,s20] = snplot(s,N,14);
title('S-N-data with estimated N(s)','FontSize',20) #title('S-N-data with estimated N(s)','FontSize',20)
set(gca,'FontSize',20) #set(gca,'FontSize',20)
if (printing==1), print -deps ../bilder/fatigue_18b.eps #if (printing==1), print -deps ../bilder/fatigue_18b.eps
end #end
wafostamp([],'(ER)') #wafostamp([],'(ER)')
disp('Block 25'),pause(pstate) #disp('Block 25'),pause(pstate)
#
#!#! Section 4.4.5 From S-N curve to fatigue life distribution ##!#! Section 4.4.5 From S-N curve to fatigue life distribution
#!#! Damage intensity as function of $\beta$ ##!#! Damage intensity as function of $\beta$
beta = 3:0.1:8; #beta = 3:0.1:8;
DRFC = cc2dam(RFC_sea,beta); #DRFC = cc2dam(RFC_sea,beta);
dRFC = DRFC/T_sea; #dRFC = DRFC/T_sea;
plot(beta,dRFC), axis([3 8 0 0.25]) #plot(beta,dRFC), axis([3 8 0 0.25])
title('Damage intensity as function of \beta') #title('Damage intensity as function of \beta')
if (printing==1), print -deps ../bilder/fatigue_19.eps #if (printing==1), print -deps ../bilder/fatigue_19.eps
end #end
wafostamp([],'(ER)') #wafostamp([],'(ER)')
disp('Block 26'),pause(pstate) #disp('Block 26'),pause(pstate)
#
#!#! Fatigue life distribution with sea load. ##!#! Fatigue life distribution with sea load.
dam0 = cc2dam(RFC_sea,beta0)/T_sea; #dam0 = cc2dam(RFC_sea,beta0)/T_sea;
[t0,F0] = ftf(e0,dam0,s20,0.5,1); #[t0,F0] = ftf(e0,dam0,s20,0.5,1);
[t1,F1] = ftf(e0,dam0,s20,0,1); #[t1,F1] = ftf(e0,dam0,s20,0,1);
[t2,F2] = ftf(e0,dam0,s20,5,1); #[t2,F2] = ftf(e0,dam0,s20,5,1);
plot(t0,F0,t1,F1,t2,F2) #plot(t0,F0,t1,F1,t2,F2)
title('Fatigue life distribution function') #title('Fatigue life distribution function')
if (printing==1), print -deps ../bilder/fatigue_20.eps #if (printing==1), print -deps ../bilder/fatigue_20.eps
end #end
wafostamp([],'(ER)') #wafostamp([],'(ER)')
disp('Block 27, last block') #disp('Block 27, last block')

291
pywafo/src/wafo/misc.py
Unescape Escape View File

@ -12,6 +12,7 @@ from numpy import (abs, amax, any, logical_and, arange, linspace, atleast_1d,atl
finfo, inf, pi, interp, isnan, isscalar, zeros, ones, finfo, inf, pi, interp, isnan, isscalar, zeros, ones,
r_, sign, unique, hstack, vstack, nonzero, where, extract) r_, sign, unique, hstack, vstack, nonzero, where, extract)
from scipy.special import gammaln from scipy.special import gammaln
from scipy.integrate import trapz, simps
import types import types
import warnings import warnings
from wafo import plotbackend from wafo import plotbackend
@ -30,7 +31,7 @@ __all__ = ['is_numlike','JITImport', 'DotDict', 'Bunch', 'printf', 'sub_dict_sel
'findoutliers', 'common_shape', 'argsreduce', 'findoutliers', 'common_shape', 'argsreduce',
'stirlerr', 'getshipchar', 'betaloge', 'gravity', 'nextpow2', 'stirlerr', 'getshipchar', 'betaloge', 'gravity', 'nextpow2',
'discretize', 'discretize2', 'pol2cart', 'cart2pol', 'meshgrid', 'ndgrid', 'discretize', 'discretize2', 'pol2cart', 'cart2pol', 'meshgrid', 'ndgrid',
'trangood', 'tranproc', 'histgrm', 'num2pistr'] 'trangood', 'tranproc', 'plot_histgrm', 'num2pistr']
def is_numlike(obj): def is_numlike(obj):
@ -59,7 +60,7 @@ class JITImport(object):
except: except:
if self._module is None: if self._module is None:
self._module = __import__(self._module_name, None, None, ['*']) self._module = __import__(self._module_name, None, None, ['*'])
assert(isinstance(self._module, types.ModuleType), 'module') #assert(isinstance(self._module, types.ModuleType), 'module')
return getattr(self._module, attr) return getattr(self._module, attr)
else: else:
raise raise
@ -1942,44 +1943,46 @@ def tranproc(x, f, x0, *xi):
if N > 4: if N > 4:
warnings.warn('Transformation of derivatives of order>4 not supported.') warnings.warn('Transformation of derivatives of order>4 not supported.')
return y #y0,y1,y2,y3,y4 return y #y0,y1,y2,y3,y4
def good_bins(data=None, range=None, num_bins=None, num_data=None, odd=False, loose=True):
def histgrm(data, n=None, odd=False, scale=False, lintype='b-'): ''' Return good bins for histogram
'''
Plot histogram
Parameters Parameters
----------- ----------
data : array-like data : array-like
the data the data
n : scalar integer range : (float, float)
approximate number of bins wanted (default depending on length(data)) minimum and maximum range of bins (default data.min(), data.max())
num_bins : scalar integer
approximate number of bins wanted (default depending on num_data=len(data))
odd : bool odd : bool
placement of bins (0 or 1) (default 0) placement of bins (0 or 1) (default 0)
scale : bool loose : bool
argument for scaling (default 0) if True add extra space to min and max
scale = 1 yields the area 1 under the histogram if False the bins are made tight to the min and max
lintype : specify color and lintype, see PLOT for possibilities.
Returns Example
------- -------
binwidth = the width of each bin >>> import wafo.misc as wm
>>> wm.good_bins(range=(0,5), num_bins=6)
Example: array([-1., 0., 1., 2., 3., 4., 5., 6.])
R=rndgumb(2,2,1,100); >>> wm.good_bins(range=(0,5), num_bins=6, loose=False)
histgrm(R,20,0,1) array([ 0., 1., 2., 3., 4., 5.])
hold on >>> wm.good_bins(range=(0,5), num_bins=6, odd=True)
x=linspace(-3,16,200); array([-1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5])
plot(x,pdfgumb(x,2,2),'r') >>> wm.good_bins(range=(0,5), num_bins=6, odd=True, loose=False)
hold off array([-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5])
''' '''
x = np.atleast_1d(data) if data:
if n is None: x = np.atleast_1d(data)
n = np.ceil(4 * np.sqrt(np.sqrt(len(x)))) num_data = len(x)
mn = x.min() mn, mx = range if range else (x.min(), x.max())
mx = x.max()
d = (mx - mn) / n * 2 if num_bins is None:
num_bins = np.ceil(4 * np.sqrt(np.sqrt(num_data)))
d = float(mx - mn) / num_bins * 2
e = np.floor(np.log(d) / np.log(10)); e = np.floor(np.log(d) / np.log(10));
m = np.floor(d / 10 ** e) m = np.floor(d / 10 ** e)
if m > 5: if m > 5:
@ -1988,11 +1991,65 @@ def histgrm(data, n=None, odd=False, scale=False, lintype='b-'):
m = 2 m = 2
d = m * 10 ** e d = m * 10 ** e
mn = (np.floor(mn / d) - 1) * d - odd * d / 2 mn = (np.floor(mn / d) - loose) * d - odd * d / 2
mx = (np.ceil(mx / d) + 1) * d + odd * d / 2 mx = (np.ceil(mx / d) + loose) * d + odd * d / 2
limits = np.arange(mn, mx, d) limits = np.arange(mn, mx+d/2, d)
return limits
def plot_histgrm(data, bins=None, range=None, normed=False, weights=None, lintype='b-'):
'''
Plot histogram
Parameters
-----------
data : array-like
the data
bins : int or sequence of scalars, optional
If an int, it defines the number of equal-width
bins in the given range (4 * sqrt(sqrt(len(data)), by default).
If a sequence, it defines the bin edges, including the
rightmost edge, allowing for non-uniform bin widths.
range : (float, float), optional
The lower and upper range of the bins. If not provided, range
is simply ``(data.min(), data.max())``. Values outside the range are
ignored.
normed : bool, optional
If False, the result will contain the number of samples in each bin.
If True, the result is the value of the probability *density* function
at the bin, normalized such that the *integral* over the range is 1.
weights : array_like, optional
An array of weights, of the same shape as `data`. Each value in `data`
only contributes its associated weight towards the bin count
(instead of 1). If `normed` is True, the weights are normalized,
so that the integral of the density over the range remains 1
lintype : specify color and lintype, see PLOT for possibilities.
Returns
-------
h : list
of plot-objects
Example
-------
>>> import pylab as plb
>>> import wafo.misc as wm
>>> import wafo.stats as ws
>>> R = ws.weibull_min.rvs(2,loc=0,scale=2, size=100)
>>> h0 = wm.plot_histgrm(R, 20, normed=True)
>>> x = linspace(-3,16,200)
>>> h1 = plb.plot(x,ws.weibull_min.pdf(x,2,0,2),'r')
See also
--------
wafo.misc.good_bins
numpy.histogram
'''
bin, limits = np.histogram(data, bins=limits, normed=scale) #, new=True) x = np.atleast_1d(data)
if bins is None:
bins = np.ceil(4 * np.sqrt(np.sqrt(len(x))))
bin, limits = np.histogram(data, bins=bins, normed=normed, weights=weights) #, new=True)
limits.shape = (-1, 1) limits.shape = (-1, 1)
xx = limits.repeat(3, axis=1) xx = limits.repeat(3, axis=1)
xx.shape = (-1,) xx.shape = (-1,)
@ -2003,9 +2060,7 @@ def histgrm(data, n=None, odd=False, scale=False, lintype='b-'):
yy[:, 0] = 0.0 # histogram yy[:, 0] = 0.0 # histogram
yy.shape = (-1,) yy.shape = (-1,)
yy = np.hstack((yy, 0.0)) yy = np.hstack((yy, 0.0))
plotbackend.plotbackend.plot(xx, yy, lintype, limits, limits * 0) return plotbackend.plotbackend.plot(xx, yy, lintype, limits, limits * 0)
binwidth = d
return binwidth
def num2pistr(x, n=3): def num2pistr(x, n=3):
''' '''
@ -2041,29 +2096,28 @@ def num2pistr(x, n=3):
xtxt = format % x xtxt = format % x
return xtxt return xtxt
def fourier(x,t,T=None,M=None,N=None, method='trapz'): def fourier(data,t=None,T=None,m=None,n=None, method='trapz'):
''' '''
Returns Fourier coefficients. Returns Fourier coefficients.
CALL: [a,b] = fourier(t,x,T,M);
Parameters Parameters
---------- ----------
x : array-like data : array-like
vector or matrix of row vectors with data points shape p x n. vector or matrix of row vectors with data points shape p x n.
t : array-like t : array-like
vector with n values indexed from 1 to N. vector with n values indexed from 1 to N.
T : real scalar T : real scalar
primitive period of signal, i.e., smallest period. (default T = t[-1]-t[0] primitive period of signal, i.e., smallest period. (default T = t[-1]-t[0]
M : scalar integer m : scalar integer
defines no of harmonics desired (default M = N) defines no of harmonics desired (default M = N)
N : scalar integer n : scalar integer
no of data points (default len(t)) no of data points (default len(t))
method : string
integration method used
Returns Returns
------- -------
a,b = Fourier coefficients size M x P a,b = Fourier coefficients size m x p
FOURIER finds the coefficients for a Fourier series representation FOURIER finds the coefficients for a Fourier series representation
of the signal x(t) (given in digital form). It is assumed the signal of the signal x(t) (given in digital form). It is assumed the signal
@ -2083,15 +2137,20 @@ def fourier(x,t,T=None,M=None,N=None, method='trapz'):
Example Example
------- -------
>>> t = linspace(0,4*T) >>> import wafo.misc as wm
>>> x = sin(t); >>> import numpy as np
>>> a, b = fourier(t, x, T=2*pi, M=5) >>> T = 2*np.pi
>>> t = np.linspace(0,4*T)
>>> x = np.sin(t)
>>> a, b = wm.fourier(x, t, T=T, m=5)
>>> (np.round(a.ravel()), np.round(b.ravel()))
(array([ 0., -0., 0., -0., 0.]), array([ 0., 4., -0., -0., 0.]))
See also See also
-------- --------
fft fft
''' '''
x = np.atleast_2d(x) x = np.atleast_2d(data)
p, n = x.shape p, n = x.shape
if t is None: if t is None:
t = np.arange(n) t = np.arange(n)
@ -2099,112 +2158,48 @@ def fourier(x,t,T=None,M=None,N=None, method='trapz'):
t = np.atleast_1d(t) t = np.atleast_1d(t)
n = len(t) if n is None else n n = len(t) if n is None else n
m = n if n is none else m m = n if n is None else m
T = t[-1]-t[0] if T is None else T T = t[-1]-t[0] if T is None else T
if method.startswith('trapz'):
intfun = trapz
# switch 0 elif method.startswith('simp'):
# case -1, intfun = simps
# % Define the vectors for computing the Fourier coefficients
# % # Define the vectors for computing the Fourier coefficients
# a = zeros(M,P); t.shape = (1,-1)
# b = zeros(M,P); a = zeros((m,p))
# a(1,:) = simpson(x); b = zeros((m,p))
# a[0] = intfun(x,t, axis=-1)
# %
# % Compute M-1 more coefficients # Compute M-1 more coefficients
# tmp = 2*pi*t(:,ones(1,P))/T; tmp = 2*pi*t/T
# % tmp = 2*pi*(0:N-1).'/(N-1); #% tmp = 2*pi*(0:N-1).'/(N-1);
# for n1 = 1:M-1, for i in range(1,m):
# n = n1+1; a[i] = intfun(x*cos(i*tmp),t, axis=-1)
# a(n,:) = simpson(x.*cos(n1*tmp)); b[i] = intfun(x*sin(i*tmp),t, axis=-1)
# b(n,:) = simpson(x.*sin(n1*tmp));
# end a = a/pi
# b = b/pi
# a = 2*a/N;
# b = 2*b/N; # Alternative: faster for large M, but gives different results than above.
# # nper = diff(t([1 end]))/T; %No of periods given
# case 0, # if nper == round(nper):
# % # N1 = n/nper
# a = zeros(M,P); # else:
# b = zeros(M,P); # N1 = n
# 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);
# #
# % # # Fourier coefficients by fft
# % Compute the dc-level (the a(0) component). # Fcof1 = 2*ifft(x(1:N1,:),[],1);
# % # Pcor = [1; exp(sqrt(-1)*(1:M-1).'*t(1))]; % correction term to get
# % 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 # % the correct integration limits
# Fcof = Fcof1(1:M,:).*Pcor(:,ones(1,P)); # Fcof = Fcof1(1:M,:).*Pcor(:,ones(1,P));
# a = real(Fcof(1:M,:)); # a = real(Fcof(1:M,:));
# b = imag(Fcof(1:M,:)); # b = imag(Fcof(1:M,:));
# end
# return return a, b

24
pywafo/src/wafo/objects.py
Unescape Escape View File

@ -100,7 +100,8 @@ class LevelCrossings(WafoData):
logcmin = logcros[icmax] logcmin = logcros[icmax]
logcros = sqrt(2 * abs(logcros - logcmin)) logcros = sqrt(2 * abs(logcros - logcmin))
logcros[0:icmax + 1] = 2 * logcros[icmax] - logcros[0:icmax + 1] logcros[0:icmax + 1] = 2 * logcros[icmax] - logcros[0:icmax + 1]
p = polyfit(self.args[10:-9], logcros[10:-9], 1) #least square fit ncr = 10
p = polyfit(self.args[ncr:-ncr], logcros[ncr:-ncr], 1) #least square fit
if self.sigma is None: if self.sigma is None:
self.sigma = 1.0 / p[0] #estimated standard deviation of x self.sigma = 1.0 / p[0] #estimated standard deviation of x
if self.mean is None: if self.mean is None:
@ -449,9 +450,10 @@ class CyclePairs(WafoData):
>>> h1 = mm.plot(marker='x') >>> h1 = mm.plot(marker='x')
''' '''
def __init__(self, *args, **kwds): def __init__(self, *args, **kwds):
self.kind = kwds.pop('kind', 'max2min') self.kind = kwds.pop('kind', 'min2max')
self.sigma = kwds.pop('sigma', None) self.sigma = kwds.pop('sigma', None)
self.mean = kwds.pop('mean', None) self.mean = kwds.pop('mean', None)
self.time = kwds.pop('time', 1)
options = dict(title=self.kind + ' cycle pairs', options = dict(title=self.kind + ' cycle pairs',
xlab='min', ylab='max', xlab='min', ylab='max',
@ -505,8 +507,8 @@ class CyclePairs(WafoData):
amp = abs(self.amplitudes()) amp = abs(self.amplitudes())
return atleast_1d([K * np.sum(amp ** betai) for betai in beta]) return atleast_1d([K * np.sum(amp ** betai) for betai in beta])
def level_crossings(self, kind='uM'): def level_crossings(self, kind='uM', intensity=False):
""" Return number of upcrossings from a cycle count. """ Return level crossing spectrum from a cycle count.
Parameters Parameters
---------- ----------
@ -516,6 +518,9 @@ class CyclePairs(WafoData):
1,'uM' : upcrossings and maxima (default). 1,'uM' : upcrossings and maxima (default).
2,'umM': upcrossings, minima, and maxima. 2,'umM': upcrossings, minima, and maxima.
3,'um' : upcrossings and minima. 3,'um' : upcrossings and minima.
intensity : bool
True if level crossing intensity spectrum
False if level crossing count spectrum
Return Return
------ ------
lc : level crossing object lc : level crossing object
@ -594,7 +599,11 @@ class CyclePairs(WafoData):
dcount = cumsum(extr[1, 0:nx]) - extr[3, 0:nx] dcount = cumsum(extr[1, 0:nx]) - extr[3, 0:nx]
elif defnr == 3: ## This are upcrossings + minima + maxima elif defnr == 3: ## This are upcrossings + minima + maxima
dcount = cumsum(extr[1, 0:nx]) + extr[2, 0:nx] dcount = cumsum(extr[1, 0:nx]) + extr[2, 0:nx]
return LevelCrossings(dcount, levels, mean=self.mean, sigma=self.sigma) ylab='Count'
if intensity:
dcount = dcount/self.time
ylab = 'Intensity [count/sec]'
return LevelCrossings(dcount, levels, mean=self.mean, sigma=self.sigma, ylab=ylab)
class TurningPoints(WafoData): class TurningPoints(WafoData):
''' '''
@ -705,7 +714,7 @@ class TurningPoints(WafoData):
""" """
if h>0: if h>0:
ind = findrfc(self.data, h, method=method) ind = findrfc(self.data, h, method=method)
data = self.data(ind) data = self.data[ind]
else: else:
data = self.data data = self.data
if data[0] > data[1]: if data[0] > data[1]:
@ -724,7 +733,8 @@ class TurningPoints(WafoData):
kind = 'max2min' kind = 'max2min'
M = data[iM:-1:2] M = 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) time = self.args[-1]-self.args[0]
return CyclePairs(M, m, kind=kind, mean=self.mean, sigma=self.sigma, time=time)
def mat2timeseries(x): def mat2timeseries(x):
""" """

Write Preview
Loading…
Cancel
Save