HERMITETR Calculate transformation, g, proposed by Winterstein Assumption: a Gaussian process, Y, is related to the non-Gaussian process, X, by Y = g(X). CALL: [g,test] = hermitetr(x,data,def); [g,test] = hermitetr(x,S,def); g = [x g(x)] a two column matrix with the transformation g(x). test = int (g(x)-x)^2 dx where int. limits is given by X. This is a measure of departure from the Gaussian model. x = a row vector with x-values. (default linspace(-5*sigma,5*sigma,501)+mean) data = [sigma skew kurt mean] is the standard deviation, skewness, kurtosis and mean of the process, respectively. skew=kurt-3=0 for a Gaussian process. This is fairly accurate if kurt>=0 and 0<=skew^2 <= 8*kurt/9 (default [1 0.16 3.04 0]) S = spectral density struct from which [sigma skew kurt mean] is calculated using spec2skew. def = 1 Winterstein et. al. (1994) parametrization (default) 2 Winterstein (1988) parametrization HERMITETR is a hermite transformation model where the transformation is chosen to be monotonic cubic polynomial, calibrated such that the first 4 moments of the transformed model G(y)=g^-1(y) match the moments of the true process. Information about the moments of the process can be obtained by site specific data, laboratory measurements or by resort to theoretical models (see spec2skew). If kurt<3 (hardening model) g(x) = xn - c3(xn^2-1) - c4*(xn^3-3*xn) where xn = (x-mean)/sigma c3 = skew/6 c4 = (kurt-3)/24 If kurt>=3 (softening model) G(y) = mean + K*sigma*[ y + c3(y^2-1) + c4*(y^3-3*y) ] where y = g(x) = G^-1(x) K = 1/sqrt(1+2*c3^2+6*c4^2) If def = 2 : c3 = skew/(6*(1+6*c4)) c4 = [sqrt(1+1.5*(kurt-3))-1]/18 If def = 1 : c3 = skew/6*(1-0.015*abs(skew)+0.3*skew^2)/(1+0.2*(kurt-3)) c4 = 0.1*((1+1.25*(kurt-3))^(1/3)-1)*c41 c41 = (1-1.43*skew^2/(kurt-3))^(1-0.1*(kurt)^0.8) NOTE: - by specifying NaN's in the data vector default values will be used. - if length(data) is shorter than the parameters needed then the default values are used for the parameters not specified. - The gaussian process in the transformed world is N(0,1) Example: Simulate a Transformed Gaussian process: Hm0=7;Tp=11; S = jonswap([],[Hm0 Tp]); g=hermitetr*Hm0/4; ys = spec2sdat(S,15000); % Simulated in the Gaussian world xs = gaus2dat(ys,g); % Transformed to the real world See also spec2skew, ochitr, lc2tr, dat2tr
Calculates the transformation by a Hermite polynomial. | |
Estimates the moments of 2'nd order non-linear waves | |
Create object or return object class. | |
Display message and abort function. | |
Linearly spaced vector. | |
Display warning message; disable or enable warning messages. |
% CHAPTER2 Modelling random loads and stochastic waves | |
% CHAPTER3 Demonstrates distributions of wave characteristics | |
% CHAPTER4 contains the commands used in Chapter 4 of the tutorial | |
Estimate transformation, g, from data. |
001 function [g ,t0]=hermitetr(x,data,def) 002 %HERMITETR Calculate transformation, g, proposed by Winterstein 003 % 004 % Assumption: a Gaussian process, Y, is related to the 005 % non-Gaussian process, X, by Y = g(X). 006 % 007 % CALL: [g,test] = hermitetr(x,data,def); 008 % [g,test] = hermitetr(x,S,def); 009 % 010 % g = [x g(x)] a two column matrix with the transformation g(x). 011 % test = int (g(x)-x)^2 dx where int. limits is given by X. This 012 % is a measure of departure from the Gaussian model. 013 % x = a row vector with x-values. 014 % (default linspace(-5*sigma,5*sigma,501)+mean) 015 % data = [sigma skew kurt mean] is the standard deviation, 016 % skewness, kurtosis and mean of the process, 017 % respectively. skew=kurt-3=0 for a Gaussian process. 018 % This is fairly accurate if kurt>=0 and 019 % 0<=skew^2 <= 8*kurt/9 (default [1 0.16 3.04 0]) 020 % S = spectral density struct from which 021 % [sigma skew kurt mean] is calculated using spec2skew. 022 % def = 1 Winterstein et. al. (1994) parametrization (default) 023 % 2 Winterstein (1988) parametrization 024 % 025 % HERMITETR is a hermite transformation model where the transformation is 026 % chosen to be monotonic cubic polynomial, calibrated such that the first 027 % 4 moments of the transformed model G(y)=g^-1(y) match the moments of 028 % the true process. Information about the moments of the process can be 029 % obtained by site specific data, laboratory measurements or by resort to 030 % theoretical models (see spec2skew). 031 % 032 % If kurt<3 (hardening model) 033 % g(x) = xn - c3(xn^2-1) - c4*(xn^3-3*xn) 034 % where 035 % xn = (x-mean)/sigma 036 % c3 = skew/6 037 % c4 = (kurt-3)/24 038 % 039 % If kurt>=3 (softening model) 040 % G(y) = mean + K*sigma*[ y + c3(y^2-1) + c4*(y^3-3*y) ] 041 % where 042 % y = g(x) = G^-1(x) 043 % K = 1/sqrt(1+2*c3^2+6*c4^2) 044 % If def = 2 : 045 % c3 = skew/(6*(1+6*c4)) 046 % c4 = [sqrt(1+1.5*(kurt-3))-1]/18 047 % If def = 1 : 048 % c3 = skew/6*(1-0.015*abs(skew)+0.3*skew^2)/(1+0.2*(kurt-3)) 049 % c4 = 0.1*((1+1.25*(kurt-3))^(1/3)-1)*c41 050 % c41 = (1-1.43*skew^2/(kurt-3))^(1-0.1*(kurt)^0.8) 051 % 052 % NOTE: - by specifying NaN's in the data vector default values will be used. 053 % - if length(data) is shorter than the parameters needed then the 054 % default values are used for the parameters not specified. 055 % - The gaussian process in the transformed world is N(0,1) 056 % 057 % Example: Simulate a Transformed Gaussian process: 058 % Hm0=7;Tp=11; 059 % S = jonswap([],[Hm0 Tp]); g=hermitetr*Hm0/4; 060 % ys = spec2sdat(S,15000); % Simulated in the Gaussian world 061 % xs = gaus2dat(ys,g); % Transformed to the real world 062 % 063 % See also spec2skew, ochitr, lc2tr, dat2tr 064 065 % References: 066 % Winterstein, S.R. (1988) 067 % 'Nonlinear vibration models for extremes and fatigue.' 068 % J. Engng. Mech., ASCE, Vol 114, No 10, pp 1772-1790 069 % 070 % Winterstein, S.R, Ude, T.C. and Kleiven, G. (1994) 071 % "Springing and slow drift responses: 072 % predicted extremes and fatigue vs. simulation" 073 % In Proc. 7th International behaviour of Offshore structures, (BOSS) 074 % Vol. 3, pp.1-15 075 076 077 078 % tested on matlab 5.2 079 % History: 080 % revised pab dec 2003 081 % revised pab 02.04.2001 082 % -Added Spectrum as a possible input 083 %revised pab 02.01.2000 084 % - added t0, check on that g is monotone 085 % - moved some code to hermitefun 086 % by pab 21.02.2000 087 088 089 data2=[1 0.16 3.04 0]; % default values 090 if nargin>=2 & ~isempty(data) 091 switch class(data) 092 case 'double', 093 if any(~isnan(data)) 094 ind=find(~isnan(data(1:min(length(data),4)))); 095 data2(ind)=data(ind); 096 end 097 case 'struct', % data is a spctral density struct 098 [skew, kurt, ma, sigma] = spec2skew(data); 099 data2 = [sigma skew, kurt, ma]; 100 otherwise 101 warning('Wrong input data') 102 end 103 end 104 sigma=data2(1); skew=data2(2); kurt=data2(3); ma=data2(4); 105 if nargin<1|isempty(x), x = linspace(-5*sigma,5*sigma,501)+ma; end 106 if nargin<3|isempty(def), def = 1;end 107 108 %skew,ga2 109 ga2 = kurt-3; 110 if ga2<0 111 c4 = ga2/24; 112 c3 = skew/6; 113 else 114 switch def 115 case 2, % Winterstein 1988 parametrization 116 if skew^2>8*(ga2+3)/9, 117 disp('warning: kurtosis too low compared to the skewness') 118 end 119 c4 = (sqrt(1+1.5*ga2)-1)/18; 120 c3 = skew/(6*(1+6*c4)); 121 otherwise, % Winterstein et. al. 1994 parametrization intended to 122 % apply for the range: 0 <= ga2 < 12 and 0<= skew^2 < 2*ga2/3 123 if skew^2>2*(ga2)/3, 124 disp('Warning: kurtosis too low compared to the skewness') 125 end 126 if (ga2 < 0)| (12 < ga2) 127 disp('Warning: kurtosis must be between 0 and 12') 128 end 129 c3 = skew/6*(1-0.015*abs(skew)+0.3*skew^2)/(1+0.2*ga2); 130 if ga2==0, 131 c4=0; 132 else 133 c41= (1-1.43*skew.^2./ga2).^(1-0.1*(ga2+3).^0.8); 134 c4 = 0.1*((1+1.25*ga2)^(1/3)-1)*c41; 135 end 136 end 137 end 138 if ~isreal(c4)|~isreal(c3) 139 error('Unable to calculate the polynomial') 140 end 141 142 if nargout>1 143 [g, t0]=hermitefun([c3 c4],x(:),ma,sigma,ga2); 144 else 145 g = hermitefun([c3 c4],x(:),ma,sigma,ga2); 146 end 147 return 148 149 150 151 152 153
Comments or corrections to the WAFO group