Makima修正Hermite插值

function [yi] = interp1TestHermite(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

% 修正Akima插值计算点的导数
len = length(x);
extendInt = 0;
slope = (y(2 : end) - y(1 : end -1))./(x(2 : end) - x(1 : end -1));

w1 = zeros(1, len-1);
w2 = zeros(1, len-1);

for i = 3 : len-2
    w1(i) = abs(slope(i+1+extendInt) - slope(i+extendInt)) + abs(slope(i+1+extendInt) + slope(i+extendInt))/2;
    w2(i) = abs(slope(i-1+extendInt) - slope(i-2+extendInt)) + abs(slope(i-1+extendInt) + slope(i-2+extendInt))/2;
end
d = zeros(1, len-1);
for i = 3 : len-2
    d(i) = w1(i)/(w1(i)+w2(i))*slope(i-1+extendInt) + w2(i)/(w1(i)+w2(i)) * slope(i+extendInt);
end

% 三次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);
    
    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;
   
%     % 已知x的间距全部为1
%     frac = xi(i)-x(xIndex);
%     fracInv = frac - 1;
%     yi(i) = ((1+2*frac)*y(xIndex) + frac * d(xIndex))*fracInv*fracInv  + ...
%             ((1-2*fracInv)*y(xIndex+1) + fracInv * d(xIndex+1))*frac*frac;
end

end

function [yi] = interp1TestHermiteMakima(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

% 修正Akima插值计算点的导数
len = length(x);
if len < 2
    error('the lenth of x can not be smaller than 2!');
end
extendInt = 2;
slope0 = (y(2 : end) - y(1 : end -1))./(x(2 : end) - x(1 : end -1));
% slope 需要向左和右各扩展两个点
slope = zeros(1, length(slope0)+2*extendInt);
slope(3 : end-2) = slope0;
if len >= 3
    slope(1) = 2*slope0(1) - slope0(2);
    slope(2) = 3*slope0(1) - 2 * slope0(2);
    slope(end-1) = 2*slope0(end) - slope0(end-1);
    slope(end) = 3*slope0(end) - 2*slope0(end-1);
end
if len == 2
    slope(1) = slope0(1);
    slope(2) = slope0(1);
    slope(end-1) = slope0(1);
    slope(end) = slope0(1);
end

w1 = zeros(1, len-1);
w2 = zeros(1, len-1);

for i = 3-extendInt : len-2+extendInt
    w1(i) = abs(slope(i+1+extendInt) - slope(i+extendInt)) + abs(slope(i+1+extendInt) + slope(i+extendInt))/2;
    w2(i) = abs(slope(i-1+extendInt) - slope(i-2+extendInt)) + abs(slope(i-1+extendInt) + slope(i-2+extendInt))/2;
end
d = zeros(1, len-1);
for i = 3-extendInt : len-2+extendInt
    d(i) = w1(i)/(w1(i)+w2(i))*slope(i-1+extendInt) + w2(i)/(w1(i)+w2(i)) * slope(i+extendInt);
end
d(isnan(d)) = 0;


% 三次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

end

clc
clear all; close all;

%% 测试只考虑中间部分的hermite插值
slen = 16;
space = 0.1;
x = 1 : slen;
y = rand(1, slen)*10;
xi = 3 : space : slen-2 ;
yi = interp1(x, y, xi, 'makima');
yis = interp1TestHermite(x, y, xi);

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

%% 测试只有两个点的makima修正的hermite插值
slen = 2;
space = 0.1;
x = 1 : slen;
y = rand(1, slen)*10;
xi = 1 : space : slen ;
yi = interp1(x, y, xi, 'makima');
yis = interp1TestHermiteMakima(x, y, xi);

figure
plot(x,y,'*-', xi, yi, 'r*-', xi, yis, 'bo-');
legend('origin data','matlab interp1-makima', 'hermite-makima(two points)');
title('interpolation hermite-makima only two point test');
mean((yis - yi))

%% 测试包含多个点的makima修正hermite插值
slen = 16;
space = 0.1;
x = 1 : slen;
y = rand(1, slen)*10;
xi = 1 : space : slen ;
yi = interp1(x, y, xi, 'makima');
yis = interp1TestHermiteMakima(x, y, xi);

figure
plot(x,y,'*-', xi, yi, 'r*-', xi, yis, 'bo-');
legend('origin data','matlab interp1-makima', 'hermite-makima');
title('interpolation hermite-makima extral four-points');
mean((yis - yi))

代码如上所示，分别测试了不外扩的makima修正三阶hermite插值、外扩4个点的makima修正三阶hermite插值(仅包含两个点，包含多个点的情况)。结果如图所示。