Home > wafo > trgauss > spec2tpdf2.m

# spec2tpdf2

## PURPOSE

Evaluates densities for various wave periods or wave lengths

## SYNOPSIS

[f] = spec2tpdf2(spec,utc,def,paramt,varargin)

## DESCRIPTION

``` SPEC2TPDF2 Evaluates densities for various wave periods or wave lengths

CALL:  F   = spec2tpdf2(S,u,def,paramt,options);

F    = density structure.
S    = spectral density structure
u    = reference level (default the most frequently crossed level).
def   = 'Tc',    gives half wave period, Tc (default).
'Tt',    gives half wave period, Tt
'Lc' and 'Lt' ditto for wave length.
paramt  = [t0 tn Nt] where t0, tn and Nt is the first value, last value
and the number of points, respectively, for which
the density will be computed. paramt= [5 5 51] implies
that the density is computed only for T=5 and
using 51 equidistant points in the interval [0,5].
options = rind-options structure containing optional parameters
controlling the performance of the integration.
See rindoptset for details.
[]    = default values are used.

SPEC2TPDF2 calculates pdf of halfperiods  Tc, Tt, Lc or Lt
in a stationary Gaussian transform process X(t),
where Y(t) = g(X(t)) (Y zero-mean Gaussian with spectrum given in S).
The transformation, g, can be estimated using LC2TR,
DAT2TR, HERMITETR or OCHITR.

Example: % The density of Tc is computed by:
S = jonswap;
paramt = [0 10 51];
f = spec2tpdf2(S,[],'Tc',paramt);
pdfplot(f)
hold on,
plot(f.x{:}, f.f+f.err,'r',f.x{:}, f.f-f.err) % estimated error bounds
hold off

## CROSS-REFERENCE INFORMATION

This function calls:
 createpdf PDF class constructor dat2gaus Transforms x using the transformation g. fwaitbar Fast display of wait bar. gaus2dat Transforms xx using the inverse of transformation g. initoptions Initializes RIND options according to speed. levels Calculates discrete levels given the parameter matrix. pdfplot Plot contents of pdf structures rind Computes multivariate normal expectations rindoptset Create or alter RIND OPTIONS structure. spec2cov2 Computes covariance function and its derivatives, alternative version spec2spec Transforms between different types of spectra wnormspec Normalize a spectral density such that m0=m2=1 clock Current date and time as date vector. close Close figure. datestr String representation of date. erf Error function. error Display message and abort function. etime Elapsed time. isequal True if arrays are numerically equal. isfield True if field is in structure array. lower Convert string to lowercase. now Current date and time as date number. num2str Convert number to string. (Fast version) pause Wait for user response. sprintf Write formatted data to string. strcmpi Compare strings ignoring case. struct2cell Convert structure array to cell array. toeplitz Toeplitz matrix.
This function is called by:

## SOURCE CODE

```001 function [f] = spec2tpdf2(spec,utc,def,paramt,varargin)
002 %SPEC2TPDF2 Evaluates densities for various wave periods or wave lengths
003 %
004 %  CALL:  F   = spec2tpdf2(S,u,def,paramt,options);
005 %
006 %        F    = density structure.
007 %        S    = spectral density structure
008 %        u    = reference level (default the most frequently crossed level).
009 %       def   = 'Tc',    gives half wave period, Tc (default).
010 %               'Tt',    gives half wave period, Tt
011 %               'Lc' and 'Lt' ditto for wave length.
012 %     paramt  = [t0 tn Nt] where t0, tn and Nt is the first value, last value
013 %               and the number of points, respectively, for which
014 %               the density will be computed. paramt= [5 5 51] implies
015 %               that the density is computed only for T=5 and
016 %               using 51 equidistant points in the interval [0,5].
017 %     options = rind-options structure containing optional parameters
018 %               controlling the performance of the integration.
019 %               See rindoptset for details.
020 %       []    = default values are used.
021 %
022 % SPEC2TPDF2 calculates pdf of halfperiods  Tc, Tt, Lc or Lt
023 % in a stationary Gaussian transform process X(t),
024 % where Y(t) = g(X(t)) (Y zero-mean Gaussian with spectrum given in S).
025 % The transformation, g, can be estimated using LC2TR,
026 % DAT2TR, HERMITETR or OCHITR.
027 %
028 % Example: % The density of Tc is computed by:
029 %    S = jonswap;
030 %    paramt = [0 10 51];
031 %    f = spec2tpdf2(S,[],'Tc',paramt);
032 %    pdfplot(f)
033 %    hold on,
034 %    plot(f.x{:}, f.f+f.err,'r',f.x{:}, f.f-f.err) % estimated error bounds
035 %    hold off
036 %
038
039 %History:
040 % Tested on : matlab 5.3
041 % revised pab 24.11.2003
042 %  updated call to rind
043 %  -removed nit and speed from input and replaced with rindoptions structure
044 % revised pab
045 % -removed the printing to screen, replaced with a call to fwaitbar
046 % -fixed some bugs: default values for speed, nit and plotflag was not
047 % properly set
048 % revised jr 16.02.2000  nit -> NIT in call of rind and f.nit = NIT;
049 % revised ir 10.02.2000 adopted to MATLAB6
050 % revised pab 23.05.2000 new name spec2tpdf2
051 % revised by Per A. Brodtkorb 20.06.1999
052 % by I. Rychlik 01.10.1998 with name wave_th1.m
053
054
055 defaultSpeed   = 9;
056 defaultoptions = initoptions(defaultSpeed);
057
058 if ((nargin==1) & (nargout <= 1) &  isequal(spec,'defaults')),
059   f = defaultoptions;
060   return
061 end
062 error(nargchk(1,inf,nargin))
063 startTime = clock;
064 if nargin<3|isempty(def)
065   def='tc';
066 end
067 if strcmpi('l',def(1))
068   spec=spec2spec(spec,'k1d');
069 elseif strcmpi('t',def(1))
070   spec=spec2spec(spec,'freq');
071 else
072   error('Unknown def')
073 end
074 switch lower(def)
075    case  {'tc','lc'},    defnr = 1; % 'tc' or 'lc'
076    case  {'tt','lt'},    defnr =-1; % 'tt' or 'lt'
077    otherwise, ,error('Unknown def')
078 end
079 [S, xl4, L0, L2, L4, L1]=wnormspec(spec);
080
081 A = sqrt(L0/L2);
082
083 if nargin<6
084   options = defaultoptions;
085 else
086   options = rindoptset(defaultoptions,varargin{:});
087 end
088
089 if isfield(spec,'tr')
090    g = spec.tr;
091 else
092    g = [];
093 end
094 if isempty(g)
095   g  = [sqrt(L0)*(-5:0.02:5)', (-5:0.02:5)'];
096 end
097 if nargin<2|isempty(utc)
098   utc_d = gaus2dat([0, 0],g); % most frequently crossed level
099   utc   = utc_d(1,2);
100 end
101
102 % transform reference level into Gaussian level
103 uu = dat2gaus([0., utc],g);
104 u  = uu(2);
105 disp(['The level u for Gaussian process = ', num2str(u)])
106
107 if nargin<4|isempty(paramt)
108   z2 = u^2/2;
109   z  = -sign(defnr)*u/sqrt(2);
110   expectedMaxPeriod = 2*ceil(2*pi*A*exp(z)*(0.5+erf(z)/2))
111   paramt = [0, expectedMaxPeriod, 51];
112 end
113
114 t0     = paramt(1);
115 tn     = paramt(2);
116 Ntime  = paramt(3);
117 t      = levels([0, tn/A, Ntime]); % normalized times
118 Nstart = max(1 + round(t0/tn*(Ntime-1)),2); % index to starting point to
119                                      % evaluate
120
121 dt    = t(2)-t(1);
122
123 nr = 2;
124 R  = spec2cov2(S,nr,Ntime-1,dt);
125
126
127 xc   = [u; u ];
128 indI = zeros(4,1);
129 Nd   = 2;
130 Nc   = 2;
131 XdInf = 100.d0*sqrt(-R(1,3));
132 XtInf = 100.d0*sqrt(R(1,1));
133
134 B_up  = [u+XtInf, XdInf, 0];
135 B_lo  = [u,    0, -XdInf];
136 INFIN = [1 1 0];
137 BIG   = zeros(Ntime+2,Ntime+2);
138 ex    = zeros(Ntime+2,1);
139 %CC    = 2*pi*sqrt(-R(1,1)/R(1,3))*exp(u^2/(2*R(1,1)));
140 %  XcScale = log(CC)
141 XcScale = log(2*pi*sqrt(-R(1,1)/R(1,3)))+(u^2/(2*R(1,1)));
142 options.xcscale = XcScale;
143
144 opt0 = struct2cell(options);
145
146 f     = createpdf;
147 f.f   = zeros(Ntime,1);
148 f.err = f.f;
149
150 h11 = fwaitbar(0,[],sprintf('Please wait ...(start at: %s)',datestr(now)));
151 for pt = Nstart:Ntime
152   Nt      = pt-Nd;
153   Ntd     = Nt+Nd;
154   indI(2) = Nt;
155   indI(3) = Nt+1;
156   indI(4) = Ntd;
157   Ntdc    = Ntd+Nc;
158   % positive wave period
159   BIG = covinput(pt,R,BIG);
160
161   [f.f(pt), f.err(pt)]= rind(BIG(1:Ntdc,1:Ntdc), ex(1:Ntdc),...
162                 B_lo,B_up,indI,xc,Nt,opt0{:});
163
164   fwaitbar(pt/Ntime,h11,sprintf('%s Ready: %d of %d',datestr(now),pt,Ntime));
165 end
166 close(h11)
167 f.err = f.err/A;
168
169
170 switch lower(def)
171  case  'tc'
172   Htxt = 'Density of Tc';
173  case  'tt'
174   Htxt = 'Density of Tt';
175  case  'lc'
176   Htxt = 'Density of Lc';
177  case  'lt'
178   Htxt = 'Density of Lt';
179 end
180
181 if strcmpi(def(1),'l')
182   xtxt = 'wave length [m]';
183 else
184   xtxt = 'period [s]';
185 end
186 Htxt = [Htxt,sprintf('_{v =%2.5g}',utc)];
187
188 f.title   = Htxt;
189 f.labx{1} = xtxt;
190 f.x{1}    = t*A;
191 f.f       = f.f/A;
192 f.options = options;
193 f.u       = utc;
194 f.elapsedTime = etime(clock,startTime);
195
196 if nargout==0
197   pdfplot(f)
198 end
199
200 return
201
202
203 function big = covinput(pt,R,big)
204 %COVINPUT
205 %
206 % CALL  big = covinput(pt,R,big)
207 %
208 % R = [R0,R1,R2] column vectors with autocovariance and its
209 %          derivatives, i.e., Ri (i=1:2) are vectors with the 1'st and 2'nd
210 %          derivatives of R0.  size Ntime x 3
211 %
212 % the order of the variables in the covariance matrix
213 % are organized as follows:
214 % For pt>1:
215 % ||X(t2)..X(ts),..X(tn-1)|| X'(t1) X'(tn)|| X(t1) X(tn) ||
216 % = [Xt                          Xd                    Xc]
217 %
218 %where
219 %
220 % Xt= time points in the indicator function
221 % Xd= derivatives
222 % Xc=variables to condition on
223
224 % Computations of all covariances follows simple rules: Cov(X(t),X(s))=r(t,s),
225 % then  Cov(X'(t),X(s))=dr(t,s)/dt.  Now for stationary X(t) we have
226 % a function r(tau) such that Cov(X(t),X(s))=r(s-t) (or r(t-s) will give the same result).
227 %
228 % Consequently  Cov(X'(t),X(s))    = -r'(s-t)    = -sign(s-t)*r'(|s-t|)
229 %               Cov(X'(t),X'(s))   = -r''(s-t)   = -r''(|s-t|)
230 %               Cov(X''(t),X'(s))  =  r'''(s-t)  =  sign(s-t)*r'''(|s-t|)
231 %               Cov(X''(t),X(s))   =  r''(s-t)   =   r''(|s-t|)
232 %               Cov(X''(t),X''(s)) =  r''''(s-t) = r''''(|s-t|)
233
234
235 %cov(Xd)
236 Sdd = -toeplitz(R([1, pt],3));
237 %cov(Xc)
238 Scc = toeplitz(R([1 pt],1));
239 %cov(Xc,Xd)
240 Scd = [0 R(pt,2); -R(pt,2) 0];
241
242
243 if pt>2,
244   %cov(Xt)
245   Stt = toeplitz(R(1:pt-2,1));  % Cov(X(tn),X(ts))  = r(ts-tn)   = r(|ts-tn|)
246   %cov(Xc,Xt)
247   Sct = R(2:pt-1,1).';          % Cov(X(tn),X(ts))  = r(ts-tn)   = r(|ts-tn|)
248   Sct = [Sct;Sct(end:-1:1)];
249   %Cov(Xd,Xt)
250   Sdt = -R(2:pt-1,2).';         % Cov(X'(t1),X(ts)) = -r'(ts-t1) = r(|s-t|)
251   Sdt = [Sdt;-Sdt(end:-1:1)];
252   N   = pt + 2;
253
254   big(1:N,1:N) = [Stt,   Sdt.',   Sct.';...
255           Sdt,   Sdd,     Scd.';...
256           Sct,   Scd,     Scc];
257  %error('covinput')
258 else
259   N = 4;
260   big(1:N,1:N) =[ Sdd,   Scd.';...
261           Scd,   Scc];
262 end
263 if  pt<0,
264   big
265   pause
266
267 end
268 return
269
270```

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