Home > wafo > trgauss > private > hermitefun.m

# hermitefun

## PURPOSE

Calculates the transformation by a Hermite polynomial.

## SYNOPSIS

[g, t0] = hermitefun(coef,x,ma,sa,gdef)

## DESCRIPTION

``` HERMITEFUN Calculates the transformation by a Hermite polynomial.
Assumption: a Gaussian process, Y, is related to the
non-Gaussian process, X, by Y = g(X).

CALL:  [g,test] = hermitefun(coef,x,ma,sa,gdef)

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 of the data from the Gaussian model.
coef   = [c3 c4], vector with polynomial coefficients, see below.
x      =  vector of x values (default linspace(-5,5,513)'*sa+ma)
ma, sa = mean and standard deviation of the process, respectively.
(default ma=0,sa=1)
gdef   = integer defining the transformation (default 1)

If gdef<0 hardening model (kurtosis < 3)
g(x) =  xn - c3(xn^2-1) - c4*(xn^3-3*xn)
where
xn = (x-ma)/sa

If gdef>=0 softening model (kurtosis >= 3)
G(y) = mean + K*sa*[ 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)

## CROSS-REFERENCE INFORMATION

This function calls:
 smooth Calculates a smoothing spline. error Display message and abort function. linspace Linearly spaced vector. num2str Convert number to string. (Fast version) polyval Evaluate polynomial. roots Find polynomial roots. sprintf Write formatted data to string. trapz Trapezoidal numerical integration. warning Display warning message; disable or enable warning messages. xor Logical EXCLUSIVE OR.
This function is called by:
 hermitetr Calculate transformation, g, proposed by Winterstein

## SOURCE CODE

```001 function [g, t0] = hermitefun(coef,x,ma,sa,gdef)
002 %HERMITEFUN Calculates the transformation by a Hermite polynomial.
003 %            Assumption: a Gaussian process, Y, is related to the
004 %                      non-Gaussian process, X, by Y = g(X).
005 %
006 % CALL:  [g,test] = hermitefun(coef,x,ma,sa,gdef)
007 %
008 %    g     = [x g(x)] a two column matrix with the transformation g(x).
009 %    test  = int (g(x)-x)^2 dx  where int. limits is given by X. This
010 %           is a measure of departure of the data from the Gaussian model.
011 %   coef   = [c3 c4], vector with polynomial coefficients, see below.
012 %   x      =  vector of x values (default linspace(-5,5,513)'*sa+ma)
013 %   ma, sa = mean and standard deviation of the process, respectively.
014 %            (default ma=0,sa=1)
015 %   gdef   = integer defining the transformation (default 1)
016 %
017 %  If gdef<0 hardening model (kurtosis < 3)
018 %     g(x) =  xn - c3(xn^2-1) - c4*(xn^3-3*xn)
019 %  where
020 %    xn = (x-ma)/sa
021 %
022 %  If gdef>=0 softening model (kurtosis >= 3)
023 %     G(y) = mean + K*sa*[ y + c3(y^2-1) + c4*(y^3-3*y) ]
024 %   where
025 %     y  = g(x) = G^-1(x)
026 %     K  = 1/sqrt(1+2*c3^2+6*c4^2)
027 %
029
030
031 % Tested on: matlab 5.3
032 % by pab
033 % - default x is now levels([-5 5 513])*sa+ma ->
034 %  better to have the discretization
035 %  represented with exact numbers, especially when calculating
036 %  derivatives of the transformation numerically.
037
038 error(nargchk(1,5,nargin))
039 if nargin<5|isempty(gdef), gdef =1 ; end
040 if nargin<4|isempty(sa), sa = 1;end
041 if nargin<3|isempty(ma), ma = 0; end
042 if nargin<2|isempty(x),  x  = linspace(-5*sa+ma,5*sa+ma,513)'; end
043
044
045 c3 = coef(1); c4 = coef(2);
046 if ~isreal(c4)|~isreal(c3)
047  error('Unable to calculate the polynomial')
048 end
049 xn = (x-ma)/sa;
050
051 if gdef<0
052   p = [-c4 -c3 1+3*c4 c3];
053   g = [x polyval(p,xn)];
054 else
055   g = [x zeros(length(xn),1)];
056   K = 1/sqrt(1+2*c3^2+6*c4^2);
057   p = K*[c4 c3 1-3*c4 -c3];%+[ 0 0 0 ma]; % polynomial coefficients
058 end
059
060
061
062 % Strip leading zeros and throw away.
063 inz= find(abs(p)>eps);
064 p = p(inz(1):end);
065 % Check if it is a strictly increasing function.
066 k  = length(p);
067 dp = [k-1:-1:1].*p(1:k-1);  % Derivative
068 r  = roots(dp);             % Find roots of the derivative
069 r(abs(imag(r))>eps)=[]; % Keep only real roots
070
071 if any(r)
072   % Check if it is possible to invert the polynomial for the given
073   % interval
074   if gdef<0
075     x00 = r;
076   else
077     x00 = sa*polyval(p,r)+ma;
078   end
079   txt1 = 'The polynomial is not a strictly increasing function.';
080   txt2 = ['The derivative of g(x) is infinite at x = ' num2str(x00')];
081   warning(sprintf('%s \n %s ',txt1,txt2))
082
083   txt2 = ['for the given interval x = ' num2str(x([1 end]).',3)];
084
085   if any(x(1)<= x00 & x00 <= x(end)),
086     cdef = 1;
087   else
088     cdef = sum( xor(x00 <= x(1) , x00 <= x(end)));
089   end
090   if mod(cdef,2),
091     txt1 = 'Unable to invert the polynomial';
092     error(sprintf('%s \n %s ',txt1,txt2))
093   end
094   txt1='However, successfully inverted the polynomial ';
095   disp(sprintf('%s \n %s ',txt1,txt2))
096 end
097
098
099 % Inverting the polynomial
100 %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
101 switch k*(gdef>=0)
102   case 0, % g(x) ready
103   case 2, % Linear
104     g(:,2) = xn;             % (x-p(2))/p(1);
105   case 3, % Quadratic: Solve c3*u^2+u-c3 = (x-ma)/(K*sa) = xn/K
106     u =  .5*(-1-sqrt(1+4*c3*(c3+xn/K)))/c3;
107     if 0,%y0>y(1)
108       g(:,2) = u/c3;             % choose the largest solution
109     else
110       g(:,2) = -(c3+xn/K)./u/c3; % choose the smallest solution
111     end
112   case 4, % Cubic: solve u + c3*(u.^2-1) + c4*(u.^3-3*u)= (x-ma)/(K*sa)
113     b = 1/(3*c4);
114     x0 = c3*b ;
115     % substitue u = z-x0  and divide by c4 => z^3 + 3*c*z+2*q0  = 0
116     c  = b-1-x0^2;
117     q0 = x0^3-1.5*b*(x0+xn/K); %
118     if any(r) & cdef ==0 % Three real roots
119       d        = sqrt(-c);
120       theta1   = acos(-q0./d^3)/3;
121       th2      = [0 -2*pi/3 2*pi/3];
122       x1       = abs(2*d*cos(theta1(ceil(length(xn)/2)) + th2)-x0);
123       [tmp ix] = min(x1);   % choose the smallest solution
124       g(:,2)   = 2*d*cos(theta1 + th2(ix))-x0;
125     else                 %Only one real root exist
126       q1 = sqrt((q0).^2+c^3);
127       % Find the real root of the monic polynomial
128       A0 = (q1-q0).^(1/3);
129       B0 = -(q1+q0).^(1/3);
130       g(:,2) = A0+B0-x0;   % real root
131       %% The other complex roots are given by
132       %x= -(A0+B0)/2+(A0-B0)*sqrt(3)/2-x0;
133       %x=-(A0+B0)/2+(A0-B0)*sqrt(-3)/2-x0;
134     end
135  case 6,
136    % old call
137    g(:,2) = xn;
138    g(:,1) = sa*polyval(p,xn)+ma;
139    %g(:,1) = sa*K*(xn + c3*(xn.^2-1) + c4*(xn.^3-3*xn))+ma;
140
141    % interpolate so that g(:,1) have equidistant spacing
142    g(:,2)=smooth(g(:,1),g(:,2),1,x,1);
143    g(:,1)=x;
144 end
145
146 if nargout>1,
147   t0 = trapz(xn,(xn-g(:,2)).^2);
148 end
149
150```

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