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)
 
  See also  hermitetr, ochitr, dat2tr

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

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 %
028 % See also  hermitetr, ochitr, dat2tr
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