LC2SDAT Simulates process with given irregularity factor and crossing spectrum CALL: L = lc2sdat(lc,N,alpha); L = two column matrix with times and values of the simulated process. lc = two column matrix with levels and crossing intensities, i.e., level-crossing spectrum.. N = number of sample points. alpha = irregularity factor, 0<alpha<1, small alpha gives irregular process.
Extracts level-crossing spectrum from data, | |
Estimate transformation, g, from observed crossing intensity. | |
Transforms process X and up to four derivatives | |
Clear current figure. | |
Error function. | |
Display message and abort function. | |
One-dimensional digital filter. | |
Hold current graph. | |
Linearly spaced vector. | |
Linear plot. | |
Graph title. | |
Trapezoidal numerical integration. |
% CHAPTER4 contains the commands used in Chapter 4 of the tutorial |
001 function process=lc2sdat(lc,N,alpha) 002 %LC2SDAT Simulates process with given irregularity factor and crossing spectrum 003 % 004 % CALL: L = lc2sdat(lc,N,alpha); 005 % 006 % L = two column matrix with times and values of the simulated 007 % process. 008 % lc = two column matrix with levels and crossing intensities, 009 % i.e., level-crossing spectrum.. 010 % N = number of sample points. 011 % alpha = irregularity factor, 0<alpha<1, small alpha gives 012 % irregular process. 013 % 014 015 % Example 016 % n = 10000; 017 % S = jonswap(7); 018 % alpha = spec2char(S,'alpha'); 019 % xs = spec2sdat(S,n); 020 % lc = dat2lc(xs); 021 % xs2 = lc2sdat(lc,n,alpha); 022 % Se = dat2spec(xs2); 023 % wspecplot(S),hold on 024 % wspecplot(Se,'r') 025 026 %History 027 % revised pab Feb2004 028 % -added example 029 % -changed order of input to conform with the other XXX2sdat routines 030 % Revised by GL, July 10, 2000: 031 % Changed call, from cross2tr to lc2tr 032 % Changed "Check the result" by removing reference to getrfc 033 034 % TODO % add a good example 035 036 error(nargchk(1,3,nargin)) 037 038 f = linspace(0,0.49999,1000); 039 rho_st = 2*sin(f*pi).^2-1; 040 tmp = alpha*asin(sqrt((1+rho_st)/2)); 041 tmp = sin(tmp).^2; 042 a2 = (tmp-rho_st)./(1-tmp); 043 y = min([a2+rho_st;1-a2]); 044 [maximum,maxidx]=max(y); 045 046 rho_st = rho_st(maxidx); 047 a2 = a2(maxidx); 048 a1 = 2*rho_st+a2-1; 049 r0 = 1; 050 r1 = -a1/(1+a2); 051 r2 = (a1^2-a2-a2^2)/(1+a2); 052 sigma2 = r0+a1*r1+a2*r2; 053 054 e = randn(N,1)*sqrt(sigma2); 055 e(1:2) = [0;0]; 056 L0 = randn(1,1); 057 L0 = [L0;r1*L0+sqrt(1-r2^2)*randn(1,1)]; 058 %Simulate the process, starting in L0 059 L = filter(1,[1 a1 a2],e,filter([1 a1 a2],1,L0)); 060 061 epsilon = 1.01; 062 min_L = min(L); 063 max_L = max(L); 064 maxi = max(abs([min_L max_L]))*epsilon; 065 mini = -maxi; 066 067 u = linspace(mini,maxi,101)'; 068 G = (1+erf(u/sqrt(2)))/2; 069 G = G.*(1-G); 070 071 x = linspace(0,r1,100)'; 072 factor1 = 1./sqrt(1-x.^2); 073 factor2 = 1./(1+x); 074 integral = zeros(size(u)); 075 for i=1:length(integral) 076 y = factor1.*exp(-u(i)*u(i).*factor2); 077 integral(i) = trapz(x,y); 078 end 079 G = G-integral/(2*pi); 080 G = G/max(G); 081 082 Z = ((u>=0)*2-1).*sqrt(-2*log(G)); 083 084 sumcr = trapz(lc(:,1),lc(:,2)); 085 lc(:,2) = lc(:,2)/sumcr; 086 mcr = trapz(lc(:,1),lc(:,1).*lc(:,2)); 087 scr = trapz(lc(:,1),lc(:,1).^2.*lc(:,2)); 088 scr = sqrt(scr-mcr^2); 089 g = lc2tr(lc,mcr,scr); 090 091 f = [u u]; 092 f(:,2) = tranproc(Z,fliplr(g)); 093 094 process = tranproc(L,f); 095 process = [(1:length(process))' process]; 096 097 %Check the result: 098 %save load.dat process -ascii 099 %[RFC TP]=getrfc; 100 %load rfc.dat 101 %param=[min(process(:,2)) max(process(:,2)) 100]; 102 %fr=cc2fr(param,RFC); 103 %NT=fr2nt(fr); 104 %mu=pickdiag(NT); 105 106 %Check the result without reference to getrfc: 107 LCe = dat2lc(process); 108 max(lc(:,2)) 109 max(LCe(:,2)) 110 111 clf 112 plot(lc(:,1),lc(:,2)/max(lc(:,2))) 113 hold on 114 plot(LCe(:,1),LCe(:,2)/max(LCe(:,2)),'-.') 115 title('Relative crossing intensity') 116 117 %% Plot made by the function funplot_4, JE 970707 118 %param = [min(process(:,2)) max(process(:,2)) 100]; 119 %plot(lc(:,1),lc(:,2)/max(lc(:,2))) 120 %hold on 121 %plot(levels(param),mu/max(mu),'--') 122 %hold off 123 %title('Crossing intensity') 124 %watstamp; 125 126 % Temporarily 127 %funplot_4(lc,param,mu); 128
Comments or corrections to the WAFO group