Home > wafo > misc > smooth.m

smooth

PURPOSE ^

Calculates a smoothing spline.

SYNOPSIS ^

[yy,coefs]= smooth(x,y,p,xx,LinExtrap,d2)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

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