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   
  
  See also  spec2cov2, wnormspec, dat2tr, dat2gaus, perioddef, wavedef

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

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 % 
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

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