样条函数拟合

样条函数拟合方法
样条函数拟合是将复杂曲线分为多段，段内用3次多项式进行拟合，同时保证分段的左右短连续且一、二次可导，以保证连接处的光滑。

Public Function FitBySpline(ByVal X() As Single, ByVal Y() As Single, ByVal SectorIndexs() As Integer) As Array ‘按照X点生成拟合后Y坐标
'SectorIndex段索引，SectorIndexs段端所指的点索引
Dim Y1(UBound(X)) As Single, i As Integer, SectorIndex As Integer, XLK As Single, Conf As Array
Conf = FitBySpline(X, Y, SectorIndexs) ‘拟合出曲线系数
For i = 0 To UBound(X)
SectorIndex = SearchSectorIndex(SearchIndex(X(i), X), SectorIndexs)
XLK = 2 * (X(i) - X(SectorIndexs(SectorIndex))) / (X(SectorIndexs(SectorIndex + 1)) - X(SectorIndexs(SectorIndex))) – 1
Y1(i) = Conf(7 * SectorIndex + 0) + Conf(7 * SectorIndex + 1) * XLK + Conf(7 * SectorIndex + 2) * (2 * XLK ^ 2 - 1) + Conf(7 * SectorIndex + 3) * (4 * XLK ^ 3 - 3 * XLK) Next
Return Y1
End Function
Private Function SearchSectorIndex(ByVal Index As Integer, ByVal SectorIndexs() As Integer) As Integer
Dim i As Integer
For i = 0 To UBound(SectorIndexs) - 1
If Index >= SectorIndexs(i) And Index <= SectorIndexs(i + 1) Then
Return i
End If
If Index >= SectorIndexs(i) And SectorIndexs(i + 1) = 0 Then '环狀
Return i
End If
Next
End Function
Private Function SearchIndex(ByVal X1 As Single, ByVal X() As Single) As Integer Dim i As Integer
For i = 0 To UBound(X)
If X1 >= X(i) And X1 <= X(i + 1) Then
Return i
End If
Next
End Function
'样条函数拟合计算，返回样条函数系数
'x--x序列
'y--序列
'SectorIndexs--分段点索引序列
'技术要求：双端在X的分段，X0和Xn是分段的首末点
Private Function FitBySpline(ByVal X() As Single, ByVal Y() As Single, ByVal SectorIndexs() As Integer) As Array
Dim SectorCount = UBound(SectorIndexs) 'M个点将差出M-1段
Dim TotalRow As Integer = (SectorCount) * 7 - 3
Dim MatrixA(TotalRow - 1, TotalRow - 1), MatrixB(TotalRow - 1) As Single
Dim CurrentSectorIndex, CurrentPointIndex, CurrentPointIndexPerSector As Integer
Dim H1, XLK As Single, PointCountPerSector As Integer
Dim A1, A2, A3, A4 As Single
For CurrentSectorIndex = 0 To SectorCount - 1
H1 = 1 / (X(SectorIndexs(CurrentSectorIndex + 1)) - X(SectorIndexs(CurrentSectorIndex)))
PointCountPerSector = SectorIndexs(CurrentSectorIndex + 1) - SectorIndexs(CurrentSectorIndex) '双端在X上，仅计入一个
For CurrentPointIndexPerSector = 0 To PointCountPerSector - 1
CurrentPointIndex = CurrentPointIndexPerSector + SectorIndexs(CurrentSectorIndex)
XLK = 2 * (X(CurrentPointIndex) - X(SectorIndexs(CurrentSectorIndex))) * H1 - 1
A1 = 1
A2 = XLK
A3 = 2 * XLK ^ 2 - 1
A4 = (4 * XLK ^ 2 - 3) * XLK
MatrixB(CurrentSectorIndex * 7 + 0) = MatrixB(CurrentSectorIndex * 7 + 0) + A1 * Y(CurrentPointIndex)
MatrixB(CurrentSectorIndex * 7 + 1) = MatrixB(CurrentSectorIndex * 7 + 1) + A2 * Y(CurrentPointIndex)
MatrixB(CurrentSectorIndex * 7 + 2) = MatrixB(CurrentSectorIndex * 7 + 2) + A3 * Y(CurrentPointIndex)
MatrixB(CurrentSectorIndex * 7 + 3) = MatrixB(CurrentSectorIndex * 7 + 3)
+ A4 * Y(CurrentPointIndex)
MatrixA(CurrentSectorIndex * 7 + 0, CurrentSectorIndex * 7 + 0) = MatrixA(CurrentSectorIndex * 7 + 0, CurrentSectorIndex * 7 + 0) + A1 ^ 2
MatrixA(CurrentSectorIndex * 7 + 1, CurrentSectorIndex * 7 + 1) = MatrixA(CurrentSectorIndex * 7 + 1, CurrentSectorIndex * 7 + 1) + A2 ^ 2
MatrixA(CurrentSectorIndex * 7 + 2, CurrentSectorIndex * 7 + 2) = MatrixA(CurrentSectorIndex * 7 + 2, CurrentSectorIndex * 7 + 2) + A3 ^ 2
MatrixA(CurrentSectorIndex * 7 + 3, CurrentSectorIndex * 7 + 3) = MatrixA(CurrentSectorIndex * 7 + 3, CurrentSectorIndex * 7 + 3) + A4 ^ 2
MatrixA(CurrentSectorIndex * 7 + 0, CurrentSectorIndex * 7 + 1) = MatrixA(CurrentSectorIndex * 7 + 0, CurrentSectorIndex * 7 + 1) + A1 * A2
MatrixA(CurrentSectorIndex * 7 + 1, CurrentSectorIndex * 7 + 2) = MatrixA(CurrentSectorIndex * 7 + 1, CurrentSectorIndex * 7 + 2) + A2 * A3
MatrixA(CurrentSectorIndex * 7 + 2, CurrentSectorIndex * 7 + 3) = MatrixA(CurrentSectorIndex * 7 + 2, CurrentSectorIndex * 7 + 3) + A3 * A4
MatrixA(CurrentSectorIndex * 7 + 0, CurrentSectorIndex * 7 + 2) = MatrixA(CurrentSectorIndex * 7 + 0, CurrentSectorIndex * 7 + 2) + A1 * A3
MatrixA(CurrentSectorIndex * 7 + 1, CurrentSectorIndex * 7 + 3) = MatrixA(CurrentSectorIndex * 7 + 1, CurrentSectorIndex * 7 + 3) + A2 * A4
MatrixA(CurrentSectorIndex * 7 + 0, CurrentSectorIndex * 7 + 3) = MatrixA(CurrentSectorIndex * 7 + 0, CurrentSectorIndex * 7 + 3) + A1 * A4
Next
MatrixA(CurrentSectorIndex * 7 + 1, CurrentSectorIndex * 7 + 0) = MatrixA(CurrentSectorIndex * 7 + 0, CurrentSectorIndex * 7 + 1)
MatrixA(CurrentSectorIndex * 7 + 2, CurrentSectorIndex * 7 + 1) = MatrixA(CurrentSectorIndex * 7 + 1, CurrentSectorIndex * 7 + 2)
MatrixA(CurrentSectorIndex * 7 + 3, CurrentSectorIndex * 7 + 2) = MatrixA(CurrentSectorIndex * 7 + 2, CurrentSectorIndex * 7 + 3)
MatrixA(CurrentSectorIndex * 7 + 2, CurrentSectorIndex * 7 + 0) = MatrixA(CurrentSectorIndex * 7 + 0, CurrentSectorIndex * 7 + 2)
MatrixA(CurrentSectorIndex * 7 + 3, CurrentSectorIndex * 7 + 1) = MatrixA(CurrentSectorIndex * 7 + 1, CurrentSectorIndex * 7 + 3)
MatrixA(CurrentSectorIndex * 7 + 3, CurrentSectorIndex * 7 + 0) = MatrixA(CurrentSectorIndex * 7 + 0, CurrentSectorIndex * 7 + 3)
If CurrentSectorIndex > 0 Then
MatrixA(CurrentSectorIndex * 7 + 0, CurrentSectorIndex * 7 - 3) = -0.5
MatrixA(CurrentSectorIndex * 7 + 1, CurrentSectorIndex * 7 - 3) = 0.5
MatrixA(CurrentSectorIndex * 7 + 1, CurrentSectorIndex * 7 - 2) = -H1
MatrixA(CurrentSectorIndex * 7 + 2, CurrentSectorIndex * 7 - 3) = -0.5
MatrixA(CurrentSectorIndex * 7 + 2, CurrentSectorIndex * 7 - 2) = 4 * H1
MatrixA(CurrentSectorIndex * 7 + 2, CurrentSectorIndex * 7 - 1) = -8 * H1 ^ 2
MatrixA(CurrentSectorIndex * 7 + 3, CurrentSectorIndex * 7 - 3) = 0.5
MatrixA(CurrentSectorIndex * 7 + 3, CurrentSectorIndex * 7 - 2) = -9 * H1
MatrixA(CurrentSectorIndex * 7 + 3, CurrentSectorIndex * 7 - 1) = 48 * H1 ^ 2
MatrixA(CurrentSectorIndex * 7 - 3, CurrentSectorIndex * 7 + 0) = -1
MatrixA(CurrentSectorIndex * 7 - 3, CurrentSectorIndex * 7 + 1) = 1
MatrixA(CurrentSectorIndex * 7 - 3, CurrentSectorIndex * 7 + 2) = -1
MatrixA(CurrentSectorIndex * 7 - 3, CurrentSectorIndex * 7 + 3) = 1
MatrixA(CurrentSectorIndex * 7 - 2, CurrentSectorIndex * 7 + 1) = -2 * H1
MatrixA(CurrentSectorIndex * 7 - 2, CurrentSectorIndex * 7 + 2) = 8 * H1
MatrixA(CurrentSectorIndex * 7 - 2, CurrentSectorIndex * 7 + 3) = -18 * H1
MatrixA(CurrentSectorIndex * 7 - 1, CurrentSectorIndex * 7 + 2) = -16 * H1 ^ 2
MatrixA(CurrentSectorIndex * 7 - 1, CurrentSectorIndex * 7 + 3) = 96 * H1 ^ 2
End If
If CurrentSectorIndex < SectorCount - 1 Then
MatrixA(CurrentSectorIndex * 7 + 0, CurrentSectorIndex * 7 + 4) = 0.5
MatrixA(CurrentSectorIndex * 7 + 1, CurrentSectorIndex * 7 + 4) = 0.5
MatrixA(CurrentSectorIndex * 7 + 1, CurrentSectorIndex * 7 + 5) = H1
MatrixA(CurrentSectorIndex * 7 + 2, CurrentSectorIndex * 7 + 4) = 0.5
MatrixA(CurrentSectorIndex * 7 + 2, CurrentSectorIndex * 7 + 5) = 4 * H1
MatrixA(CurrentSectorIndex * 7 + 2, CurrentSectorIndex * 7 + 6) = 8 * H1 ^ 2
MatrixA(CurrentSectorIndex * 7 + 3, CurrentSectorIndex * 7 + 4) = 0.5
MatrixA(CurrentSectorIndex * 7 + 3, CurrentSectorIndex * 7 + 5) = 9 * H1
MatrixA(CurrentSectorIndex * 7 + 3, CurrentSectorIndex * 7 + 6) = 48 * H1 ^ 2
MatrixA(CurrentSectorIndex * 7 + 4, CurrentSectorIndex * 7 + 0) = 1
MatrixA(CurrentSectorIndex * 7 + 4, CurrentSectorIndex * 7 + 1) = 1
MatrixA(CurrentSectorIndex * 7 + 4, CurrentSectorIndex * 7 + 2) = 1
MatrixA(CurrentSectorIndex * 7 + 4, CurrentSectorIndex * 7 + 3) = 1
MatrixA(CurrentSectorIndex * 7 + 5, CurrentSectorIndex * 7 + 1) = 2 * H1
MatrixA(CurrentSectorIndex * 7 + 5, CurrentSectorIndex * 7 + 2) = 8 * H1
MatrixA(CurrentSectorIndex * 7 + 5, CurrentSectorIndex * 7 + 3) = 18 * H1
MatrixA(CurrentSectorIndex * 7 + 6, CurrentSectorIndex * 7 + 2) = 16 * H1 ^ 2
MatrixA(CurrentSectorIndex * 7 + 6, CurrentSectorIndex * 7 + 3) = 96 * H1 ^ 2
End If
Next
Dim CGsSolveLine As New CGsSolveLine()
If CGsSolveLine.LEGauss(TotalRow, MatrixA, MatrixB) Then '全选主元消去法求解线性方程组
Return MatrixB
End If
End Function
Public Class CGsSolveLine
'n-－数组维数
'DblA－－输入二维数组
'DblB－－输入一维数组，返回结果
'全选主元消去法求解线性方程组
'Note:A,B数组同时为工作数组，遭到破坏
Public Function LEGauss(ByVal n As Integer, ByVal DblA(,) As Single, ByVal DblB() As Single) As Boolean
Dim i, j, k As Integer
Dim nIs As Integer, nJs(n - 1) As Integer
Dim d, t As Single
'开始求解
For k = 0 To n - 2
d = 0
'归一
For i = k To n - 1
For j = k To n - 1
t = Math.Abs(DblA(i, j))
If t > d Then
d = t
nJs(k) = j
nIs = i
End If
Next
Next
'无解，返回
If d + 1 = 1 Then
LEGauss = False
Exit Function
End If
'消元
If nJs(k) <> k Then
For i = 0 To n - 1
t = DblA(i, k)
DblA(i, k) = DblA(i, nJs(k))
DblA(i, nJs(k)) = t
Next
End If
If nIs <> k Then
For j = k To n - 1
t = DblA(k, j)
DblA(k, j) = DblA(nIs, j)
DblA(nIs, j) = t
Next
t = DblB(k)
DblB(k) = DblB(nIs)
DblB(nIs) = t
End If
d = DblA(k, k)
For j = k + 1 To n - 1
DblA(k, j) = DblA(k, j) / d
Next
DblB(k) = DblB(k) / d
For i = k + 1 To n - 1
For j = k + 1 To n - 1
DblA(i, j) = DblA(i, j) - DblA(i, k) * DblA(k, j) Next
DblB(i) = DblB(i) - DblA(i, k) * DblB(k)
Next
Next
d = DblA(n - 1, n - 1)
'无解，返回
If Math.Abs(d) + 1 = 1 Then
Return False
End If
'回代
DblB(n - 1) = DblB(n - 1) / d
For i = n - 2 To 0 Step -1
t = 0
For j = i + 1 To n - 1
t = t + DblA(i, j) * DblB(j)
Next
DblB(i) = DblB(i) - t
Next
'调整解得顺序
nJs(n - 1) = n - 1
For k = n - 1 To 0 Step -1
If nJs(k) <> k Then
t = DblB(k)
DblB(k) = DblB(nJs(k))
DblB(nJs(k)) = t
End If
Next
Return True
End Function
End Class。