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
Calculates a smoothing spline. | |
Display message and abort function. | |
Linearly spaced vector. | |
Convert number to string. (Fast version) | |
Evaluate polynomial. | |
Find polynomial roots. | |
Write formatted data to string. | |
Trapezoidal numerical integration. | |
Display warning message; disable or enable warning messages. | |
Logical EXCLUSIVE OR. |
Calculate transformation, g, proposed by Winterstein |
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
Comments or corrections to the WAFO group