SPEC2CMAT Joint intensity matrix for cycles (max,min)-, rainflow- and (crest,trough) CALL: f = spec2cmat(S,u,def,paramt,paramu,nit); f = pdf (density structure) of crests (trough) heights S = spectral density structure u = reference level (default the most frequently crossed level). def = 'Mm' : gives maximum and the following minimum height. 'rfc' : gives maximum and the rainflow minimum height. 'AcAt': gives (crest,trough) heights (this option needs more work). paramt = [0 tn Nt] defines discretization of half period: tn is the longest period considered while Nt is the number of points, i.e. (Nt-1)/tn is the sampling frequnecy. paramt=[0 10 51] implies that the halfperiods are considered at 51 linearly spaced points in the interval [0,10], i.e. sampling frequency is 5 Hz. paramu = [u v N] defines discretization of maxima and minima ranges: u is the lowest minimum considered, v the heighest maximum and N is the number of levles (u,v) included. nit = 0,...,9. Dimension of numerical integration (only positive nit are allowed). (default nit=1). [] = default values are used. The model for loads is a stationary Gaussian transformed process X(t), where Y(t) = g(X(t)) is a zero-mean Gaussian with spectrum, S. Note: algorithm uses Markov Chain approximation to the sequence of turning points in Y. Example: % The intensity matrix of rainflow cycles is computed by: S = jonswap; L0 = spec2mom(S,1); paramu = [sqrt(L0)*[-4 4] 41]; frfc = spec2cmat(S,[],'rfc',[],paramu); See also spec2mmtpdf spec2cov, wnormspec, dat2tr, dat2gaus, wavedef
PDF class constructor | |
Transforms x using the transformation g. | |
Transforms xx using the inverse of transformation g. | |
Calculates discrete levels given the parameter matrix. | |
Rainflow matrix given a Markov matrix of a Markov chain of turning points | |
Frequencies of upcrossing troughs and crests using Markov chain of turning points. | |
Calculates quantile levels which encloses P% of PDF | |
Transforms process X and up to four derivatives | |
Calculates joint density of minimum and following maximum | |
Normalize a spectral density such that m0=m2=1 | |
Current date and time as date vector. | |
Display message and abort function. | |
Elapsed time. | |
True if field is in structure array. | |
Convert string to lowercase. | |
Convert number to string. (Fast version) | |
Display warning message; disable or enable warning messages. |
% CHAPTER1 demonstrates some applications of WAFO | |
% CHAPTER4 contains the commands used in Chapter 4 of the tutorial | |
Script to computer exercises 3 |
001 function [f, fmm] = spec2cmat(spec,utc,def,paramt,paramu,nit) 002 %SPEC2CMAT Joint intensity matrix for cycles (max,min)-, rainflow- and (crest,trough) 003 % 004 % CALL: f = spec2cmat(S,u,def,paramt,paramu,nit); 005 % 006 % f = pdf (density structure) of crests (trough) heights 007 % S = spectral density structure 008 % u = reference level (default the most frequently crossed level). 009 % def = 'Mm' : gives maximum and the following minimum height. 010 % 'rfc' : gives maximum and the rainflow minimum height. 011 % 'AcAt': gives (crest,trough) heights (this option needs 012 % more work). 013 % paramt = [0 tn Nt] defines discretization of half period: tn is 014 % the longest period considered while Nt is the number of 015 % points, i.e. (Nt-1)/tn is the sampling frequnecy. 016 % paramt=[0 10 51] implies that the halfperiods are 017 % considered at 51 linearly spaced points in the interval 018 % [0,10], i.e. sampling frequency is 5 Hz. 019 % paramu = [u v N] defines discretization of maxima and minima ranges: 020 % u is the lowest minimum considered, v the heighest 021 % maximum and N is the number of levles (u,v) included. 022 % nit = 0,...,9. Dimension of numerical integration (only 023 % positive nit are allowed). (default nit=1). 024 % [] = default values are used. 025 % 026 % 027 % The model for loads is a stationary Gaussian transformed process X(t), 028 % where Y(t) = g(X(t)) is a zero-mean Gaussian with spectrum, S. 029 % 030 % Note: algorithm uses Markov Chain approximation to the sequence of 031 % turning points in Y. 032 % 033 % Example: % The intensity matrix of rainflow cycles is computed by: 034 % S = jonswap; 035 % L0 = spec2mom(S,1); 036 % paramu = [sqrt(L0)*[-4 4] 41]; 037 % frfc = spec2cmat(S,[],'rfc',[],paramu); 038 % 039 % See also spec2mmtpdf spec2cov, wnormspec, dat2tr, dat2gaus, wavedef 040 041 042 % Tested on : matlab 5.3 043 % History: by I. Rychlik 01.10.1998 with name minmax.m 044 % bounds by I.R. 02.01.2000. 045 % Revised by es 000322. Made call with directional spectrum possible. 046 % revised by ir 000612. Help and plots improved. 047 % revised by IR removing error in transformation 29 VI 2000 048 % revised by I.R. 01.20.2001 Change name minmax to wminmax 049 % revised by I.R. 6 II 2001 adapted for MATLAB 6 050 % revised pab 30nov2003 051 052 % TODO % AcAt option needs more work 053 startTime = clock; 054 [S, xl4, L0, L2, L4, L1]=wnormspec(spec); 055 056 057 A = sqrt(L0/L2); 058 SCIS=0; 059 if nargin<6|isempty(nit) 060 nit=1; 061 elseif nit<0 062 warning('Only postive nit allowed') 063 nit=1; 064 end 065 066 if isfield(spec,'tr') 067 g = spec.tr; 068 else 069 g = []; 070 end 071 if isempty(g) 072 g = [sqrt(L0)*(-5:0.02:5)', (-5:0.02:5)']; 073 end 074 S.tr = g; 075 076 if nargin<7|isempty(speed) 077 speed=5; 078 end 079 if nargin<2|isempty(utc) 080 utc_d=gaus2dat([0, 0],g); % most frequent crossed level 081 utc=utc_d(1,2); 082 end 083 084 % transform reference level into Gaussian level 085 uu = dat2gaus([0., utc],g); 086 u = uu(2); 087 disp(['The level u for Gaussian process = ', num2str(u)]) 088 089 090 if nargin<4|isempty(paramt) 091 distanceBetweenExtremes = 5*pi*sqrt(L2/L4); %(2.5 * mean distance between extremes) 092 paramt = [0 distanceBetweenExtremes,41]; 093 end 094 t0 = paramt(1); 095 tn = paramt(2); 096 Ntime = paramt(3); 097 t = levels([0, tn/A, Ntime]); % normalized times 098 099 100 if nargin<3|isempty(def) 101 defnr=0; 102 else 103 switch lower(def) 104 case 'mm', defnr = 1; 105 case 'rfc', defnr = 0; 106 case 'acat', defnr =-1; 107 otherwise, error('Unknown def') 108 end 109 end 110 111 112 nr = 4; 113 114 if nargin<5|isempty(paramu) 115 paramu=[-4*sqrt(L0), 4*sqrt(L0), 41]; 116 end 117 %Transform amplitudes to Gaussian levels: 118 h = levels(paramu); 119 h = reshape(h,length(h),1); 120 Nx = length(h); 121 der = ones(Nx,1); 122 hg = tranproc([h, der],g); 123 der = abs(hg(:,2)); 124 hg = hg(:,1); % Gaussian level 125 126 %if exist('h.in'), delete h.in, end 127 %fid=fopen('h.in','wt'); 128 %fprintf(fid,'%12.10f\n',hg); 129 %fclose(fid); 130 131 132 133 %paru=paramu; 134 % paru(1:2)=paru(1:2)/sqrt(L0); 135 fmM = wminmax(S,nit,paramu,t); 136 137 138 Htxt=['Joint density of crest and trough']; 139 labx='crest [m]'; 140 laby='trough [m]'; 141 if (defnr==0) 142 Htxt = ['Joint density of maximum and rainflow minimum']; 143 labx='max [m]'; 144 laby='rainflow min [m]'; 145 end 146 f=createpdf; 147 f.title=Htxt; 148 f.nit=nit; 149 f.speed=speed; 150 f.SCIS=SCIS; 151 f.labx{1}=labx; 152 f.labx{2}=laby; 153 f.x{1}=h; 154 fmm=createpdf; 155 fmm.title = ['Joint density of maximum and minimum']; 156 fmm.labx{1}='max [m]'; 157 fmm.labx{2}='min [m]'; 158 fmm.nit=nit; 159 fmm.speed=speed; 160 fmm.SCIS=SCIS; 161 dh=h(2)-h(1); 162 f.x{1}=h; 163 f.x{2}=h; 164 fmm.x{1}=h; 165 fmm.x{2}=h; 166 ftmp=fmM; 167 168 for i=1:Nx 169 ftmp(:,i)=dh*dh*ftmp(:,i);%.*der*der(i);%* sqrt(-R(1,6)/R(1,4))/2/pi; 170 end 171 fmm.f=ftmp; 172 173 if (defnr==0) 174 f.f=fliplr(mctp2rfc(fliplr(ftmp)));%* sqrt(-R(1,6)/R(1,4))/2/pi; 175 fmm.f=ftmp; 176 end 177 if (defnr==1) 178 f.f=ftmp; 179 end 180 if (defnr==-1) 181 f.f = fliplr(mctp2tc(fliplr(ftmp),utc,paramu)); 182 index1 = find(f.x{1}>0); 183 index2 = find(f.x{2}<0); 184 f.f = flipud(f.f(index2,index1)); 185 f.x{1} = f.x{1}(index1); 186 f.x{2} = abs(flipud(f.x{2}(index2))); 187 end 188 [f.cl,f.pl] = qlevels(f.f,[10, 30, 50, 70, 90, 95, 99, 99.9],f.x{1},f.x{2}); 189 190 f.elapsedTime = etime(clock,startTime); 191 192
Comments or corrections to the WAFO group