polyharmonic升维函数matlab代码

Polyharmonic升维函数在matlab中的代码如下：

function [y] = polyharmonic(x, p, d)
% POLYHARMONIC Polyharmonic function for up-sampling
%
%   y = polyharmonic(x,p) computes values of the polyharmonic function
%   of order p evaluated at points x, where x is a column vector and y is
%   returned as a column vector. This function can be used to up-sample
%   signals by inserting zeros between samples and then applying a
%   polyharmonic interpolant.
%
%   y = polyharmonic(x,p,d) is similar but returns the pth derivative
%   of the polyharmonic function.
%
%   Examples:
%       % Up-sample a signal using polyharmonic interpolation
%       x = (1:8)';
%       y = sin(x);
%       xi = linspace(1,8,33)';
%       yi = polyharmonic(xi,2);
%       plot(x,y,'o',xi,yi)
%
%       % Compute the second derivative of the polyharmonic function
%       x = (1:8)';
%       y = sin(x);
%       yp = polyharmonic(x,2,1);
%       ypp = polyharmonic(x,2,2);
%       plot(x,y,x,yp,x,ypp)
%
%   See also INTERP1, SPLINE.

% Author:      Peter J. Acklam
% Time-stamp:  2004-09-28 09:49:35 +0200
% E-mail:      pjacklam@online.no
% URL:         http://home.online.no/~pjacklam

if nargin < 2
   error('Too few input arguments.');
end

if nargin < 3
   d = 0;
end

if ~isvector(x)
   error('X must be a vector.');
end

if p ~= round(p) || p < 1
   error('P must be a positive integer.');
end

if d ~= round(d) || d < 0 || d > p
   error('D must be an integer between 0 and P.');
end

x = x(:);
n = length(x);
y = zeros(n,p);

% Compute distance matrix.
dx = x(:,ones(1,n)) - x(:,ones(1,n)).';
if d == 0
   y = dx.^p .* log(dx + eps);
else
   switch d
      case 1, y = p * dx.^(p-2);
      case 2, y = p * (p-2) * dx.^(p-4);
      case 3, y = p * (p-2) * (p-4) * dx.^(p-6);
      case 4, y = p * (p-2) * (p-4) * (p-6) * dx.^(p-8);
      otherwise
         q = 1;
         for k = p-d+1:p
            q = q * k;
         end
         y = q * dx.^(p-d);
   end
end

% Compute weights.
w = y \ eye(n);

% Compute function values.
y = w(dx.^p);

% Return appropriate derivative.
if d > 0
   for k = 1:d
      y = y + (-1)^k * nchoosek(d,k) * polyharmonic(x,p-d+k-1);
   end
end

if nargout == 0
   plot(x,y)
   clear y
end