ppcreate.m

算荧光衰减时间或者寿命的软件.matlab语言编写.

function fh=ppcreate(x,y,pptype,P)
%PPCREATE Create 1-D Piecewise Polynomial Function.
% PPfun=PPCREATE(X,Y) creates a 1-D spline piecewise polynomial from the
% data in X and Y, and returns a function handle PPfun for evaluating the
% piecewise polynomial.
% PPfun=PPCREATE(X,Y,Type,P) creates a 1-D piecewise polynomial of a
% specified Type from the data in X and Y, and returns a function handle
% PPfun for evaluating the piecewise polynomial. P is an array of
% additional parameters required for some piecewise polynomial types.
%
% Type is a character string defining the type of piecewise polynomial to
% be created:
%
% Type        Description
% 'pchip'     same as MATLAB PCHIP function, no P required.
% 'spline'    same as MATLAB SPLINE function, no P required.
% 'notaknot'  same as MATLAB SPLINE function, no P required.
% 'extrap'    same as MATLAB SPLINE function, no P required.
% 'natural'   spline, y''=0 at end points, no P required.
% 'parabolic' spline, first and last polys are parabolic, no P required.
% 'periodic'  spline, y' and y'' match at end points, no P required.
% 'aperiodic' spline, y' opposite signs, y'' equal at end points, no P.
% 'clamped'   spline, P is two element vector of slopes y' at end points.
% 'curvature' spline, P is two element vector of y'' at the end points.

% PPfun(X) evaluates the piecewise polynomial at the values in X.
% PPfun('pp') returns the piecewise polynomial structure contained in PPfun.
%
% See also: SPLINE, PCHIP, PPHELP

% D.C. Hanselman, University of Maine, Orono, ME 04469
% [e-mail]
% Mastering MATLAB 7
% 2005-10-19

if nargin<2
   error('At Least Two Inputs are Required.')
end
if ~isreal(x)
   error('X Must Contain Real Values.')
end
[x,idx]=sort(x(:)); % put x in increasing order
nx=length(x);
if nx~=numel(y);
   error('X and Y Must Contain the Same Number of Elements.')
end
y=reshape(y(idx),size(x));
if nx<4
   error('At Least 4 Data Points are Required.')
end
if nargin==2
   pptype='spline';
end
if ~ischar(pptype)
   error('Type Must be a Character String.')
end
if strncmpi(pptype,'c',1)
   if nargin<4
      error('Parameter Vector P Required for Chosen Piecewise Polynomial.')
   end
   if ~isnumeric(P) || numel(P)~=2 || ~isreal(P)
      error('Parameter P Must be a Real Two Element Numeric Vector.')
   end
end
   
n=nx-1;    % number of splines
H=diff(x); % differences in x
if any(H==0)
   error('X Values Must be Distinct.')
end

if strncmpi(pptype,'pc',2)   % MATLAB pchip
   pp=pchip(x,y);
else                         % splines
   
   D=diff(y)./H;                   % slopes
   V=[0;6*diff(D);0];              % right side vector
   B=[0;2*(H(1:n-1)+H(2:n));0];    % diagonal elements
   A=[0;H(2:n)];                   % subdiagonal elements
   C=A;                            % superdiagonal elements

   switch lower(pptype(1:2))
   case {'sp' 'no' 'ex'}  % MATLAB spline
      B(1)=H(2);
      B(2)=B(2)+H(1)+H(1)*H(1)/H(2);
      B(n)=B(n)+H(n)+H(n)*H(n)/H(n-1);
      B(nx)=H(n-1);
      C(1)=-H(1)-H(2);
      C(2)=C(2)-H(1)*H(1)/H(2);
      C(n)=0;
      A(n)=-H(n-1)-H(n);
      A(n-1)=A(n-1)-H(n)*H(n)/H(n-1);
      V(1)=0;
      V(nx)=0;   
   case 'na' % natural spline
      B(1)=1;
      B(nx)=1;
      C(n)=0;
      A(n)=0;
      V(1)=0;
      V(nx)=0;
   case 'pa' % parabolic spline
      B(1)=1;
      B(2)=B(2)+H(1);
      B(n)=B(n)+H(n);
      B(nx)=-1;
      C(1)=-1;
      C(n)=0;
      A(n)=1;
      V(1)=0;
      V(nx)=0;
   case 'pe' % periodic spline
      B(1)=1;
      B(nx)=H(n)/H(1);
      C(1)=1/2;
      A(1)=H(1);
      A(n)=B(nx)/2;
      V(1)=3*(D(1)-D(n))/H(1);
      V(nx)=V(1);
   case 'ap' % aperiodic spline
      B(1)=1;
      B(nx)=-H(n)/H(1);
      C(1)=1/2;
      A(1)=H(1);
      A(n)=B(nx)/2;
      V(1)=3*(D(1)+D(n))/H(1);
      V(nx)=V(1);
   case 'cl' % clamped spline
      B(1)=1;
      B(2)=B(2)-H(1)/2;
      B(n)=B(n)-H(n)/2;
      B(nx)=1;
      C(1)=0.5;
      C(n)=0;
      A(n)=0.5;
      V(1)=3*(D(1)-P(1))/H(1);
      V(2)=V(2)-3*(D(1)-P(1));
      V(n)=V(n)-3*(P(2)-D(n));
      V(nx)=3*(P(2)-D(n))/H(n);
   case 'cu' % curvature specified spline
      B(1)=1;
      B(nx)=1;
      C(n)=0;
      A(n)=0;
      V(1)=P(1);
      V(2)=V(2)-H(1)*P(1);
      V(n)=V(n)-H(n)*P(2);
      V(nx)=P(2);
   otherwise
         error('Unknown Type Requested: %s',pptype)
   end
   i=[2:nx 1:nx  1:n ];                             % row indices for A;B;C
   j=[1:n  1:nx  2:nx];                          % column indices for A;B;C
   ABC=sparse(i,j,[A;B;C],nx,nx,3*(nx+1));           % create sparse matrix

   switch lower(pptype(1:2))
   case {'sp' 'no' 'ex'}                           % poke in extra elements
      ABC(1,3)=H(1);
      ABC(nx,n-1)=H(n);
   case 'pe'
      ABC(1,nx)=H(n)/H(1);
      ABC(1,n)=ABC(1,nx)/2;
      ABC(nx,1)=1;
      ABC(nx,2)=1/2;
   case 'ap'
      ABC(1,nx)=-H(n)/H(1);
      ABC(1,n)=ABC(1,nx)/2;
      ABC(nx,1)=1;
      ABC(nx,2)=1/2;
   end
   sflag=spparms('autommd');
   spparms('autommd',0);   % no reordering required
   m=ABC\V;                % find solution
   spparms('autommd',sflag);

   s(:,4)=y(1:n);                    % x^3 coefficients
   s(:,3)=D-H.*(2*m(1:n)+m(2:nx))/6; % x^2 coefficients
   s(:,2)=m(1:n)/2;                  % x^1 coefficients
   s(:,1)=diff(m)./(6*H);            % x^0 coefficients

   pp.form = 'pp'; % create pp structure
   pp.breaks = x.';
   pp.coefs = s;
   pp.pieces = n;
   pp.order = 4;
   pp.dim = 1;
end
fh=@(x) ppeval(pp,x);

%--------------------------------------------------------------------------
function yi=ppeval(pp,xi)
if ischar(xi) && strcmpi(xi,'pp')
   yi=pp;
else
   yi=ppval1(pp,xi);
end