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 See also spec2cov2, wnormspec, dat2tr, dat2gaus, perioddef, wavedef
PDF class constructor | |
Transforms x using the transformation g. | |
Fast display of wait bar. | |
Transforms xx using the inverse of transformation g. | |
Initializes RIND options according to speed. | |
Calculates discrete levels given the parameter matrix. | |
Plot contents of pdf structures | |
Computes multivariate normal expectations | |
Create or alter RIND OPTIONS structure. | |
Computes covariance function and its derivatives, alternative version | |
Transforms between different types of spectra | |
Normalize a spectral density such that m0=m2=1 | |
Current date and time as date vector. | |
Close figure. | |
String representation of date. | |
Error function. | |
Display message and abort function. | |
Elapsed time. | |
True if arrays are numerically equal. | |
True if field is in structure array. | |
Convert string to lowercase. | |
Current date and time as date number. | |
Convert number to string. (Fast version) | |
Wait for user response. | |
Write formatted data to string. | |
Compare strings ignoring case. | |
Convert structure array to cell array. | |
Toeplitz matrix. |
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 % 037 % See also spec2cov2, wnormspec, dat2tr, dat2gaus, perioddef, wavedef 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
Comments or corrections to the WAFO group