Home > wafo > wsim > lc2sdat.m

# lc2sdat

## PURPOSE

Simulates process with given irregularity factor and crossing spectrum

## SYNOPSIS

process=lc2sdat(lc,N,alpha)

## DESCRIPTION

``` 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.```

## CROSS-REFERENCE INFORMATION

This function calls:
 dat2lc Extracts level-crossing spectrum from data, lc2tr Estimate transformation, g, from observed crossing intensity. tranproc Transforms process X and up to four derivatives clf Clear current figure. erf Error function. error Display message and abort function. filter One-dimensional digital filter. hold Hold current graph. linspace Linearly spaced vector. plot Linear plot. title Graph title. trapz Trapezoidal numerical integration.
This function is called by:
 Chapter4 % CHAPTER4 contains the commands used in Chapter 4 of the tutorial

## SOURCE CODE

```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
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:
099 %[RFC TP]=getrfc;
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```

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group

Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003