pchip-三阶Hermite插值

使用MATLAB测试pchip方法，其中Function pchipslopes是从MATLAB处拷贝而来。

function [yi] = interp1TestHermitePchip(x, y, xi)
% 测试makima插值
% 输入 - x: 距离, y: 函数值 , xi 待插值点
% 输出 - yi: 待插值点的函数值
% 不考虑越界的情况

% 保证一维输入且x与y的长度一致
sizex = (size(x));
sizey = (size(y));
if ~(sizex(1)*sizex(2)==sizex(1) || sizex(1)*sizex(2)==sizex(2)) || ~(sizey(1)*sizey(2)==sizey(1) || sizey(1)*sizey(2)==sizey(2))
    error('the input x or y must be one dimension!')
end

if length(x) ~= length(y)
    error('the lenth of x or y must be same!');
end

% 插值计算点的导数
len = length(x);
if len < 2
    error('the lenth of x can not be smaller than 2!');
end

slope = (y(2 : end) - y(1 : end -1))./(x(2 : end) - x(1 : end -1));

d = pchipslopes(x,y,slope);


% 三次hermite插值
yi = zeros(1, length(xi));
for i = 1 : length(xi)
    xIndex = 0;
    for j = 1 : 1 : length(x)
        if xi(i) - x(j) >= 0
            xIndex = j;
        end
    end
    fprintf(2, 'xIndex = %d\n', xIndex);
    
    if xIndex < 1
        xIndex = 1;
    end
    if xIndex >= len
        xIndex = len - 1;
    end
    
    yi(i) = ((1+2*(xi(i)-x(xIndex))/(x(xIndex+1)-x(xIndex)))*y(xIndex) + (xi(i)-x(xIndex))*d(xIndex))*((xi(i)-x(xIndex+1))/(x(xIndex)-x(xIndex+1)))^2+...
       ((1+2*(xi(i)-x(xIndex+1))/(x(xIndex)-x(xIndex+1)))*y(xIndex+1) + (xi(i)-x(xIndex+1))*d(xIndex+1))*((xi(i)-x(xIndex))/(x(xIndex+1)-x(xIndex)))^2;
    
end


function d = pchipslopes(x,y,del)
% copy from MATLAB
% d = pchipslopes(x,y,del) computes the first derivatives, d(k) = P'(x(k)).

%  Special case n=2, use linear interpolation.

   n = length(x);
   if n==2  
      d = repmat(del(1),size(y));
      return
   end

%  Slopes at interior points.
%  d(k) = weighted average of del(k-1) and del(k) when they have the same sign.
%  d(k) = 0 when del(k-1) and del(k) have opposites signs or either is zero.

   d = zeros(size(y));
  
   k = find(sign(del(1:n-2)).*sign(del(2:n-1)) > 0);
  
   h = diff(x);
   hs = h(k)+h(k+1);
   w1 = (h(k)+hs)./(3*hs);
   w2 = (hs+h(k+1))./(3*hs);
   dmax = max(abs(del(k)), abs(del(k+1)));
   dmin = min(abs(del(k)), abs(del(k+1)));
   d(k+1) = dmin./conj(w1.*(del(k)./dmax) + w2.*(del(k+1)./dmax));

%  Slopes at end points.
%  Set d(1) and d(n) via non-centered, shape-preserving three-point formulae.

   d(1) = ((2*h(1)+h(2))*del(1) - h(1)*del(2))/(h(1)+h(2));
   if sign(d(1)) ~= sign(del(1))
      d(1) = 0;
   elseif (sign(del(1)) ~= sign(del(2))) && (abs(d(1)) > abs(3*del(1)))
      d(1) = 3*del(1);
   end
   d(n) = ((2*h(n-1)+h(n-2))*del(n-1) - h(n-1)*del(n-2))/(h(n-1)+h(n-2));
   if sign(d(n)) ~= sign(del(n-1))
      d(n) = 0;
   elseif (sign(del(n-1)) ~= sign(del(n-2))) && (abs(d(n)) > abs(3*del(n-1)))
      d(n) = 3*del(n-1);
   end

% testPchip.m文件
clc
% clear all; close all;

%% 测试6个点的三次样条插值
slen = 12;
space = 0.1;
x = 1 : slen;
y = rand(1, slen)*10;
xi = 0 : space : slen+1;
yi = pchip(x, y, xi);

yis = interp1TestHermitePchip(x, y, xi);

yiPchip = interp1(x,y,xi,'pchip');

mean(yi - yis)
mean(yiPchip-yi)

figure
plot(x,y,'*-', xi, yi, 'r*-', xi, yis, 'bo-');
legend('origin data','matlab pchip','test pchip');
title('interpolation')

显示的图像为：