Home > wafo > trgauss > hermitetr.m

# hermitetr

## PURPOSE

Calculate transformation, g, proposed by Winterstein

## SYNOPSIS

[g ,t0]=hermitetr(x,data,def)

## DESCRIPTION

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

## CROSS-REFERENCE INFORMATION

This function calls:
 hermitefun Calculates the transformation by a Hermite polynomial. spec2skew Estimates the moments of 2'nd order non-linear waves class Create object or return object class. error Display message and abort function. linspace Linearly spaced vector. warning Display warning message; disable or enable warning messages.
This function is called by:
 Chapter2 % CHAPTER2 Modelling random loads and stochastic waves Chapter3 % CHAPTER3 Demonstrates distributions of wave characteristics Chapter4 % CHAPTER4 contains the commands used in Chapter 4 of the tutorial dat2tr Estimate transformation, g, from data.

## SOURCE CODE

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

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