Akima 插值拟合算法VBA版

本文参考silent_dusbin的akima 插值拟合算法 Python/C++/C版本  改写为VBA版，共需要的同学使用。

Attribute VB_Name = "Alima"
Option Explicit
Option Base 0
Function AKimaResult(NewX As Double, Xrange As Range, Yrange As Range) As Double
    Dim x() As Double, y() As Double
    Dim n As Integer, i
    n = Xrange.Cells.Count
    ReDim x(n - 1), y(n - 1)
    For i = 0 To n - 1
        x(i) = Xrange.Cells(i + 1, 1)
        y(i) = Yrange.Cells(i + 1, 1)
    Next
    AKimaResult = AkimaInterpolation(NewX, x, y)
End Function

Function AkimaInterpolation(NewX As Double, xData() As Double, yData() As Double) As Double
    Dim s(0 To 3) As Double 's[0]~s[3]分别是0~3次项的项系数
    Dim Index As Integer, dX As Integer 'Index是索引值
    Dim u(0 To 4) As Double '斜率
    Dim Weights(0 To 3) As Double '斜率变化权重关系,u[0]~u[1]->weights[0],u[1]~u[2]->weights[1],以此类推
    Dim p As Double '左斜率变化甲醛平均
    Dim q As Double '右斜率变化加权平均
    Dim t As Double 't表示区间内的与最近插值点的距离，相比传统的Akima，没有做归一化
    
    Dim n As Integer, i As Integer
    
    n = UBound(xData) - LBound(xData) + 1
    If (UBound(yData) - LBound(yData) + 1) <> n Then
        MsgBox "给定x、y数据数量不等，本函数无法拟合"
        Exit Function
    End If
    '少于5个点不给拟合
    If n < 5 Then
        MsgBox "给定数据少于5个，本函数不能拟合"
        AkimaInterpolation = 0
        Exit Function
    End If
    
    '重新处理X值
    'xData(0) = 0
    'For i = 1 To n - 1
    '    xData(i) = xData(i) - xData(0)
    'Next
    'NewX = NewX - xData(0)
    'xData(0) = 0
    
    For i = 0 To 3
        s(i) = 0
    Next
    
    '寻找当前输入NewX最近的xData位置
    If NewX < xData(1) Then
        Index = 0
    ElseIf (NewX >= xData(n - 1)) Then
        Index = n - 2
    Else
        '二分法查找
        Dim iLeft As Integer, iRight As Integer, iMid As Integer
        iLeft = 1
        iRight = n
        iMid = 1
        While ((iLeft - iRight) <> 1) And ((iLeft - iRight) <> -1) And (iLeft <> iRight)
            iMid = (iLeft + iRight) / 2
            If NewX < xData(iMid - 1) Then
                iRight = iMid
            Else
                iLeft = iMid
            End If
        Wend
        Index = iLeft - 1
        '下面是遍历查找
        'For i = 0 To n - 1
        '    If (NewX >= xData(i)) And (NewX <= xData(i + 1)) Then
        '        Index = i
        '        Exit For
        '    End If
        'Next
    End If
    
    If Index >= n - 1 Then Index = n - 2
    
    '计算斜率-------------------------------------------------
    u(2) = (yData(Index + 1) - yData(Index)) / (xData(Index + 1) - xData(Index)) '当前斜率
    
    If Index <= 1 Then
        u(3) = (yData(Index + 2) - yData(Index + 1)) / (xData(Index + 2) - xData(Index + 1))
        If Index = 1 Then
            u(1) = (yData(1) - yData(0)) / (xData(1) - xData(0))
            u(0) = 2# * u(1) - u(2)
            If n = 4 Then
                u(4) = 2# * u(3) - u(2)
            Else
                u(4) = (yData(4) - yData(3)) / (xData(4) - xData(3))
            End If
        Else
            u(1) = 2# * u(2) - u(3)
            u(0) = 2# * u(1) - u(2)
            u(4) = (yData(3) - yData(2)) / (xData(3) - xData(2))
        End If
    ElseIf Index >= n - 3 Then
        u(1) = (yData(Index) - yData(Index - 1)) / (xData(Index) - xData(Index - 1))
        If Index = n - 3 Then
            u(3) = (yData(n - 1) - yData(n - 2)) / (xData(n - 1) - xData(n - 2))
            u(4) = 2# * u(3) - u(2)
            If n = 4 Then
                u(0) = 2# * u(1) - u(2)
            Else
                u(0) = (yData(Index - 1) - yData(Index - 2)) / (xData(Index - 1) - xData(Index - 2))
            End If
        Else
            u(3) = 2# * u(2) - u(1)
            u(4) = 2# * u(3) - u(2)
            u(0) = (yData(Index - 1) - yData(Index - 2)) / (xData(Index - 1) - xData(Index - 2))
        End If
    Else
        u(0) = (yData(Index - 1) - yData(Index - 2)) / (xData(Index - 1) - xData(Index - 2)) '左左斜率
        u(1) = (yData(Index) - yData(Index - 1)) / (xData(Index) - xData(Index - 1))         '左侧斜率
        u(3) = (yData(Index + 2) - yData(Index + 1)) / (xData(Index + 2) - xData(Index + 1)) '右侧斜率
        u(3) = (yData(Index + 3) - yData(Index + 2)) / (xData(Index + 3) - xData(Index + 2)) '右右斜率
    End If
    
    '计算权重--------------------------------------------------------------
    Weights(0) = Abs(u(3) - u(2)) '右侧斜率变化权重
    Weights(1) = Abs(u(1) - u(0)) '左左侧斜率变化权重
    '下面主要判断浮点数是否等于0
    If (Weights(0) + 1# = 1#) And (Weights(1) + 1# = 1#) Then 'if abs(s(0))<0.00000001 and abs(s(1))<0.00000001
        p = (u(1) + u(2)) / 2 '直接求权重平均
    Else
        p = (Weights(0) * u(1) + Weights(1) * u(2)) / (Weights(0) + Weights(1))
        '展开u[1]*(w[0]/w[0]+w[1])+u[2]*(w[1]/w[0]+w[1])相当于在配置p中斜率u[1]和u[2]的占比
    End If
    
    Weights(2) = Abs(u(4) - u(3)) '右右侧斜率变化权重
    Weights(3) = Abs(u(2) - u(1)) '左侧斜率变化权重
    If (Weights(2) + 1# = 1#) And (Weights(3) + 1# = 1#) Then
        q = (u(2) + u(3)) / 2
    Else
        q = (Weights(2) * u(2) + Weights(3) * u(3)) / (Weights(2) + Weights(3))
    End If
    '//展开u[2]*(w[2]/w[2]+w[3])+u[3]*(w[3]/w[2]+w[3])相当于在配置q中斜率u[2]和u[3]的占比
    '//从上面可以看出在配置p和q斜率加权平均的时候，u[0]和u[4]不直接影响p q,
    '//而是通过weight[2] weight[1]影响u[1] u[2] u[3]间接影响p q,u[2]是当前斜率，u[1] u[3]是左右临侧斜率

    '//s[0]~s[3]分别是0次项到3次项系数,为什么是这些系数如果有读者找到相关文献请在评论区告诉我
    s(0) = yData(Index)
    s(1) = p                                    '这里可选p或q
    dX = xData(Index + 1) - xData(Index)        '就是微分中的dx
    s(2) = (3 * u(2) - 2 * p - q) / dX          '如果s(1)选择了q，这里的p和q要调换
    s(3) = (q + p - 2 * u(2)) / (dX * dX)       '如果s(1)选择了q，这里的p和q要调换
    t = NewX - xData(Index)
    AkimaInterpolation = s(0) + s(1) * t + s(2) * t * t + s(3) * t * t * t
End Function

 源代码下载：Akima样条插值.bas