SMOOTH Calculates a smoothing spline. CALL: yy = smooth(x,y,p,xx,def,v) x = x-coordinates of data. (vector) y = y-coordinates of data. (vector) p = [0...1] is a smoothing parameter: 0 -> LS-straight line 1 -> cubic spline interpolant xx = the x-coordinates in which to calculate the smoothed function. yy = the calculated y-coordinates of the smoothed function if xx is given. If xx is not given then yy contain the pp-form of the spline, for use with PPVAL. def = 0 regular smoothing spline (default) 1 a smoothing spline with a constraint on the ends to ensure linear extrapolation outside the range of the data v = variance of each y(i) (default ones(length(X),1)) Given the approximate values y(i) = g(x(i))+e(i) of some smooth function, g, where e(i) is the error. SMOOTH tries to recover g from y by constructing a function, f, which minimizes p * sum (Y(i) - f(X(i)))^2/d2(i) + (1-p) * int (f'')^2 The call pp = smooth(x,y,p) gives the pp-form of the spline, for use with PPVAL. Example:% x = linspace(0,1).'; y = exp(x)+1e-1*randn(size(x)); pp = smooth(x,y,.9); plot(x,y,x,smooth(x,y,.99,x),'g',x,ppval(pp,x),'k',x,exp(x),'r') See also lc2tr, dat2tr, ppval
Difference and approximate derivative. | |
Display message and abort function. | |
Make piecewise polynomial. | |
Evaluate piecewise polynomial. | |
Sparse matrix formed from diagonals. |
Estimate transformation, g, from observed CDF. | |
Estimate transformation, g, from data. | |
Smooths the conditional DIST2D distribution parameters. | |
Smooths the conditional DIST2D distribution parameters. | |
Calculates the transformation by a Hermite polynomial. | |
Joint (Scf,Hd) PDF for nonlinear waves with a JONSWAP spectra. | |
Joint (Scf,Hd) PDF for linear waves with JONSWAP spectra. | |
Joint (Vcf,Hd) PDF for linear waves with a JONSWAP spectrum. | |
Joint (Vcf,Hd) PDF for linear waves with a JONSWAP spectrum. | |
Estimate transformation, g, from observed crossing intensity. | |
Estimate transformation, g, from observed crossing intensity, version2. | |
Wave height, Hd, distribution parameters for Ochi-Hubble spectra. | |
Joint (Scf,Hd) PDF for linear waves with Ochi-Hubble spectra. | |
Joint (Scf,Hd) PDF linear waves in space with Ochi-Hubble spectra. | |
Joint (Vcf,Hd) PDF for lineare waves with Ochi-Hubble spectra. | |
Transform joint T-H density to V-H density | |
Joint (Scf,Hd) PDF for nonlinear waves with Torsethaugen spectra. | |
Joint (Scf,Hd) PDF for linear waves with Torsethaugen spectra. | |
Joint (Scf,Hd) PDF for linear waves with Torsethaugen spectra. | |
Joint (Scf,Hd) PDF for linear waves in space with Torsethaugen spectra. | |
Joint (Vcf,Hd) PDF for linear waves with Torsethaugen spectra. | |
Joint distribution (pdf) of crest front velocity and wave height: |
001 function [yy,coefs]= smooth(x,y,p,xx,LinExtrap,d2) 002 %SMOOTH Calculates a smoothing spline. 003 % 004 % CALL: yy = smooth(x,y,p,xx,def,v) 005 % 006 % x = x-coordinates of data. (vector) 007 % y = y-coordinates of data. (vector) 008 % p = [0...1] is a smoothing parameter: 009 % 0 -> LS-straight line 010 % 1 -> cubic spline interpolant 011 % xx = the x-coordinates in which to calculate the smoothed function. 012 % yy = the calculated y-coordinates of the smoothed function if 013 % xx is given. If xx is not given then yy contain the 014 % pp-form of the spline, for use with PPVAL. 015 % def = 0 regular smoothing spline (default) 016 % 1 a smoothing spline with a constraint on the ends to 017 % ensure linear extrapolation outside the range of the data 018 % v = variance of each y(i) (default ones(length(X),1)) 019 % 020 % Given the approximate values 021 % 022 % y(i) = g(x(i))+e(i) 023 % 024 % of some smooth function, g, where e(i) is the error. SMOOTH tries to 025 % recover g from y by constructing a function, f, which minimizes 026 % 027 % p * sum (Y(i) - f(X(i)))^2/d2(i) + (1-p) * int (f'')^2 028 % 029 % The call pp = smooth(x,y,p) gives the pp-form of the spline, 030 % for use with PPVAL. 031 % 032 % Example:% 033 % x = linspace(0,1).'; 034 % y = exp(x)+1e-1*randn(size(x)); 035 % pp = smooth(x,y,.9); 036 % plot(x,y,x,smooth(x,y,.99,x),'g',x,ppval(pp,x),'k',x,exp(x),'r') 037 % 038 % See also lc2tr, dat2tr, ppval 039 040 041 %References: 042 % Carl de Boor (1978) 043 % 'Practical Guide to Splines' 044 % Springer Verlag 045 % Uses EqXIV.6--9, pp 239 046 047 % tested on: Matlab 4.x 5.x 048 %History: 049 % revised pab 26.11.2000 050 % - added example 051 % - fixed a bug: n=2 is now possible 052 % revised by pab 21.09.99 053 % - added d2 054 % revised by pab 29.08.99 055 % -new extrapolation: ensuring that 056 % the smoothed function has contionous derivatives 057 % in the first and last knot 058 % modified by Per A. Brodtkorb 23.09.98 059 % secret option forcing linear extrapolation outside the ends when p>0 060 % used in lc2tr 061 062 if (nargin<5)|(isempty(LinExtrap)), 063 LinExtrap=0; %do not force linear extrapolation in the ends (default) 064 end 065 066 067 [xi,ind]=sort(x(:)); 068 n = length(xi); 069 y = y(:); 070 if n<2, 071 error('There must be >=2 data points.') 072 elseif any(diff(xi)<=0), 073 error('Two consecutive values in x can not be equal.') 074 elseif n~=length(y), 075 error('x and y must have the same length.') 076 end 077 078 if nargin<6|isempty(d2), 079 d2 = ones(n,1); %not implemented yet 080 else 081 d2=d2(:); 082 end 083 084 085 yi=y(ind); %yi=yi(:); 086 087 dx = diff(xi); 088 dydx = diff(yi)./dx; 089 090 if (n==2) % straight line 091 coefs=[dydx yi(1)]; 092 else 093 if LinExtrap & n==3, p = 0;end % Force LS-fit 094 dx1=1./dx; 095 Q = spdiags([dx1(1:n-2) -(dx1(1:n-2)+dx1(2:n-1)) dx1(2:n-1)],0:-1:-2,n,n-2); 096 D = spdiags(d2,0,n,n); % The variance 097 R = spdiags([dx(1:n-2) 2*(dx(1:n-2)+dx(2:n-1)) dx(2:n-1)],-1:1,n-2,n-2); 098 099 QQ = (6*(1-p))*(Q.'*D*Q)+p*R; 100 % Make sure Matlab uses symmetric matrix solver 101 u = 2*((QQ+QQ')\diff(dydx)); % faster than u=QQ\(Q'*yi); 102 ai = yi-6*(1-p)*D*diff([0;diff([0;u;0]).*dx1;0]); % faster than yi-6*(1-p)*Q*u 103 104 % The piecewise polynominals are written as 105 % fi=ai+bi*(x-xi)+ci*(x-xi)^2+di*(x-xi)^3 106 % where the derivatives in the knots according to Carl de Boor are: 107 % ddfi = 6*p*[0;u] = 2*ci; 108 % dddfi = 2*diff([ci;0])./dx = 6*di; 109 % dfi = diff(ai)./dx-(ci+di.*dx).*dx = bi; 110 111 ci = 3*p*[0;u]; 112 113 114 if LinExtrap & p~=0& n>3, %Forcing linear extrapolation in the ends 115 ci([2, end]) = 0; 116 % New call 117 % fixing the coefficients so that we have continous 118 % derivatives everywhere 119 ai(1) = -(ai(3)-ai(2))*dx(1)/dx(2) +ai(2)+ ci(3)*dx(1)*dx(2)/3; 120 ai(n) = (ai(n-1)-ai(n-2))*dx(n-1)/dx(n-2) +ai(n-1)+ ci(n-2)*dx(n-2)*dx(n-1)/3; 121 end 122 123 di = diff([ci;0]).*dx1/3; 124 bi = diff(ai).*dx1-(ci+di.*dx).*dx; 125 coefs = [di ci bi ai(1:n-1)]; 126 end 127 128 pp = mkpp(xi,coefs); 129 if (nargin<4)|(isempty(xx)), 130 yy = pp; 131 else 132 yy = ppval(pp,xx); 133 end 134 135 136 return 137 % old call: linear extrapolation 138 % crude the derivatives are not necessarily continous 139 % where we force ci and di to zero 140 if LinExtrap, %forcing linear extrapolation in the ends 141 ci([2, end])=0; % could be done better 142 di([1 end])=0; 143 end 144 % Old call for LinExtrap 145 if 0 %LinExtrap & n>6 & (p~=0), %forcing linear extrapolation in the ends 146 Q = Q(2:end-1,2:end-1); 147 D = D(2:end-1,2:end-1); 148 R = R(2:end-1,2:end-1); 149 % new call : so ci(1:2)=ci(n-1:n)=di(1)=di(n)=0 150 151 QQ=(6*(1-p))*(Q.'*D*Q)+p*R; 152 % Make sure Matlab uses symmetric matrix solver 153 u=2*((QQ+QQ')\diff(dydx(2:n-2))); %faster than 2*(QQ+QQ.'')\(Q'*yi); 154 155 % The piecewise polynominals are written as 156 % fi=ai+bi*(x-xi)+ci*(x-xi)^2+di*(x-xi)^3 157 ai=yi(2:n-1)-(6*(1-p)*D*diff([0;diff([0;u;0])./dx(2:n-2);0]));%faster than yi-6*(1-p)*D*Q*u ; 158 ci=3*p*[0;0;u;0]; 159 % fixing the coefficients so that we have continous 160 % derivatives everywhere 161 a1=-(ai(2)-ai(1))*dx(1)/dx(2) +ai(1)+ ci(3)*dx(1)*dx(2)/3; 162 an=(ai(n-2)-ai(n-3))*dx(n-1)/dx(n-2) +ai(n-2)+ ci(n-2)*dx(n-2)*dx(n-1)/3; 163 ai=[a1;ai; an]; 164 165 else 166 167 end 168 169 170
Comments or corrections to the WAFO group