spaps.m

演示matlab曲线拟和与插直的基本方法

function [sp,values,rho] = spaps(x,y,tol,varargin)
%SPAPS Smoothing spline.
%
%   [SP,VALUES] = SPAPS(X,Y,TOL)  returns the B-form and, if asked, the
%   values at X, of a cubic smoothing spline  f  to the data X, Y.
%   The smoothing spline approximates, at the data site X(j), the given
%   data value Y(:,j), j=1:length(X). The data values may be scalars,
%   vectors, matrices, or even ND-arrays.
%   Data points with the same site are replaced by their (weighted) average,
%   and the tolerance TOL is reduced accordingly.
%
%   The smoothing spline  f  minimizes the roughness measure
%
%       F(D^M f) := integral |D^M f(t)|^2 dt  on  min(X) < t < max(X)
%
%   over all functions  f  for which the error measure
%
%       E(f) :=  sum_j { W(j)*|Y(:,j) - f(X(j))|^2 : j=1,...,length(X) }
%
%   is no bigger than TOL.
%   
%   Here,  X  and  Y  are the result of replacing any data points with the
%   same site by their average, with its weight the sum of the corresponding
%   weights, and the given TOL is reduced correspondingly.
%   D^M f  denotes the M-th derivative of  f .
%   The weights W are chosen so that  E(f)  is the Composite Trapezoid Rule
%   approximation for  F(y-f).
%   The integer M is chosen to be 2, hence  D^M f  is the second derivative
%   of  f . For this choice of M,  f  is a CUBIC smoothing spline.
%   This is so because  f  is constructed as the unique minimizer of
%
%                     rho*E(f) + F(D^2 f) ,
%
%   with the smoothing parameter RHO so chosen that  E(f)  equals TOL.
%   Hence, FN2FM(SP,'pp') should be (up to roundoff) the same as the output
%   from CPAPS(X,Y,RHO/(1+RHO)).
%
%   [SP,VALUES,RHO] = SPAPS(X,Y,...) also returns RHO.
%
%   When TOL is 0, the variational interpolating cubic spline is obtained.
%   If TOL is negative, then  -TOL  is taken as the value of RHO to be used.
%
%   When the data values Y(:,j) are not just scalars but of size  d , then,
%   for each  t , also  f(t)  and  D^M f(t)  are of size  d , hence have
%   altogether  prod(d)  entries, in which case the expressions
%   |Y(:,j) - f(X(j))|^2 and  |D^M f(t)|^2  denote the sum of squares of
%   these  prod(d)  entries.
%
%   For functions of several variables, the data sites must be gridded,
%   as discussed below. For scattered data, try TPAPS.
%
%   Example:
%
%      x = linspace(0,2*pi,21); y = sin(x) + (rand(1,21)-.5)*.2;
%      sp = spaps(x,y, (.05)^2*(x(end)-x(1)) );
%
%   returns a fit to the noisy data that is expected to be quite close to
%   the underlying noisefree data since the latter come from a slowly
%   varying function and since the TOL used is of the size appropriate for
%   the size of the noise.
%
%   SPAPS(X,Y,TOL,W)
%   SPAPS(X,Y,TOL,M)
%   SPAPS(X,Y,TOL,W,M)
%   SPAPS(X,Y,TOL,M,W)
%   all let you provide W and/or M explicitly, with  M  presently restricted
%   to  M = 1 (linear) and M = 3 (quintic) (in addition to the default, M = 2).
%
%   You can also vary the smoothness requirement by having TOL be a sequence
%   rather than a scalar. In that case, the roughness measure is taken to be
%
%                    integral lambda(t) (D^M f(t))^2 dt ,
%
%   with  lambda  the piecewise constant function with breaks  X  whose
%   value on the interval  (X(i-1) .. X(i))  is TOL(i), all i, while TOL(1)
%   continues to be taken as the required tolerance, TOL.
%
%   Example:
%
%      sp1 = spaps(x,y, [(.025)^2*(x(end)-x(1)),ones(1,10),repmat(.1,1,10)] );
%
%   uses the same data and tolerance as the earlier example, but chooses the
%   roughness weight to be only .1 in the right half of the interval and
%   gives, correspondingly, a rougher but better fit there.
%   A plot showing both examples for comparison could now be obtained by
%
%      fnplt(sp); hold on, fnplt(sp1,'r'), plot(x,y,'ok'), hold off
%      title('Two cubic smoothing splines')
%      xlabel('the red one has reduced smoothness requirement in right half.')
%
%   SPAPS({X1,...,Xr},Y,...) provides a smoothing spline to data values
%   Y  on the r-dimensional rectangular grid specified by the  r  vectors
%   X1, ..., Xr , and these vectors may be of different lengths.
%   Now  Y  has size  [d,length(X1), ..., length(Xr)] .  If  Y  is only an
%   r-dimensional array, but of size [length(X1), ..., length(Xr)], then
%   d  is taken to be [1], i.e., the function is scalar-valued.
%   In this case of gridded data, the optional argument M is either an
%   integer or else an r-vector of integers (from the set  {1,2,3} ), and
%   the optional argument W is expected to be a cell array of length  r
%   with W{i} either empty (to get the default choice) or else a vector of
%   the same length as X{i}, i=1:r.
%   Also, if, in this case, TOL is not a cell-array, then it must either
%   be an r-vector specifying the r required tolerances, or else its first
%   entry is used as the required tolerance in each of the r variables.
%
%   Example:
%      x = -2:.2:2; y=-1:.25:1; [xx,yy] = ndgrid(x,y);
%      z = exp(-(xx.^2+yy.^2)) + (rand(size(xx))-.5)/30;
%      sp = spaps({x,y},z,8/(60^2));  fnplt(sp)
%
%   produces the picture of a smooth approximant to noisy data from a
%   bivariate function.
%   Use of MESHGRID here instead of NDGRID would produce an error.
%
%   See also SPAPI, SPAP2, CSAPS, TPAPS.

%   Copyright 1987-2005 C. de Boor and The MathWorks, Inc.
%   $Revision: 1.30.4.2 $

w = []; m = [];
for j=4:nargin
   arg = varargin{j-3};
   if ~isempty(arg)
      if iscell(arg)
         if length(arg{1})==1, m = cat(2,arg{:});
         else                  w = arg;
         end
      else
         if length(arg)==1, m = arg;
         else               w = arg;
         end
      end
   end
end
if isempty(m), m = 2; end

if iscell(x)          % we are dealing with multivariate, gridded data

   r = length(x); % can't use m here since it's used as the smoothing order
   sizey = size(y);
   if length(sizey)<r
     error('SPLINES:SPAPS:wrongsizeY', ...
          ['If X is a cell-array of length m, then Y must have', ...
            ' at least m dimensions.']), end
   if length(sizey)==r, % grid values of a scalar-valued function
     if issparse(y), y = full(y); end
     sizey = [1 sizey]; y = reshape(y, sizey);
   end

   sizeval = sizey(1:end-r); sizey = [prod(sizeval), sizey(end-r+(1:r))];
   y = reshape(y, sizey);

   if ~iscell(tol)
      tol = tol(:).';
      if length(tol)~=r, tol = repmat(tol(1),1,r); end
      tol = num2cell(tol);
   end
   if length(m)==1,  m = repmat(m,1,r); end
   if isempty(w), w = cell(1,r); end

   v = y; sizev = sizey;
   for i=r:-1:1   % carry out coordinatewise smoothing
      if nargout>2
        [sp,ignored,rho{i}] =  ...
                   spaps1(x{i}, reshape(v,prod(sizev(1:r)),sizev(r+1)),...
                          tol{i}, w{i},m(i));
        [t,a,n] = spbrk(sp);
      else
        [t,a,n] = spbrk(spaps1(x{i}, reshape(v,prod(sizev(1:r)),sizev(r+1)),...
                          tol{i}, w{i},m(i)));
      end
      knots{i} = t; sizev(r+1) = n; v = reshape(a,sizev);
      if r>1
         v = permute(v,[1,r+1,2:r]); sizev(2:r+1) = sizev([r+1,2:r]);
      end
   end
   % At this point, V contains the tensor-product B-spline coefficients,
   % and KNOTS contains the various knot sequences.
   % It remains to package the information:
   sp = spmak(knots, v); % v cannot have trailing singleton dims, hence no need
                         % to use the optional third input argument.
   if length(sizeval)>1, sp = fnchg(sp,'dz',sizeval); end

   if nargout>1, values = fnval(sp,x); end

else             % we have univariate data
   if nargout>1
      [sp,values,rho] = spaps1(x,y,tol,w,m);
   else
      sp = spaps1(x,y,tol,w,m);
   end
end

function [sp,values,rho] = spaps1(x,y,tol,w,m)
%SPAPS1 Univariate smoothing spline.

[xi,yi,sizeval,w,origint,tol,tolred] = chckxywp(x,y,max(2,m),w,tol,'adjtol');

if tol(1)>=0 % we may have to reduce the tolerance, because of data points
             % with the same site but different values ...
   tol(1) = tol(1)-tolred;
   if tol(1)<0
      warning('SPLINES:SPAPS:toltoolow', ...
           ['The specified tolerance cannot be met due to the fact that\n', ...
	   ' some data points have the same data site.']);
      tol(1) = 0;
   end
end   

dx = diff(xi); n = size(xi,1); xi = reshape(xi,1,n);
yd = size(yi,2);

if n==m % the smoothing spline is the interpolating polynomial
   rho = 0;
   values = yi.'; sp = spmak(augknt(xi([1 n]),m),yi(ones(m,1),:).');
elseif tol(1)==0||tol(1)==(-inf)  % return the interpolating variational spline;
   rho = inf; values = yi.';
   switch m
   case 1
      sp = spmak(augknt(xi,2),yi.');
   case 2
      sp = fn2fm(csape(xi,yi.','variational'),'B-');
   case 3
      % here's a quick fix on getting the quintic variational interpolating