Introduction
This is my attempt at producing an Excel function to interpolate data using splines. A good introduction to spline interpolation is, once again, given by wikipedia.
The algorithm is taken from the book 'Numerical Analysis, 7th Edition, Burden and Faires'. The macro seems to give the same answers as the example in the book so I'm fairly confident that it mostly works.
The macro takes:
- The x coordinates of a set of points
- The y coordinates of a set of points
- The x coordinates of the points that you want to interpolate
- An optional parameter of the gradient at point zero
- An optional parameter of the gradient at point n
If the optional parameters are missing then free interpolation is performed, otherwise a clamped interpolation is done. The macro outputs the interpolated Y coordinates.
Also you will need to create a class module called Point2D. The code for this class is shown below.
The Code
Public Function SplineInterpolation(dataRangeX As range, dataRangeY As range, interpX As range, Optional gradient0 As Variant, Optional gradientN As Variant) As Variant
'Spline interpolation - algorithm from Numerical Analysis 7th Edition, Burden & Faires
Dim Data() As Point2D
Dim h() As Double
Dim alpha() As Double
Dim l() As Double
Dim mu() As Double
Dim z() As Double
Dim i As Integer
Dim j As Integer
Dim numPoints As Integer
Dim c() As Double
Dim d() As Double
Dim b() As Double
Dim OutputRows As Long
Dim OutputCols As Long
Dim output() As Double
Dim vertData As Boolean
Dim vertOutput As Boolean
Dim vertInterp As Boolean
Dim x1 As Double
Dim y1 As Double
Dim fpo As Double
Dim fpn As Double
Dim jStop As Integer
'Set up output variables
With Application.Caller
OutputRows = .Rows.Count
OutputCols = .Columns.Count
End With
ReDim output(1 To OutputRows, 1 To OutputCols)
'Is the output range verticle or horizontal
If OutputRows > OutputCols Then
vertOutput = True
Else
vertOutput = False
End If
'Is the data range verticle or horizontal
If dataRangeX.Rows.Count > dataRangeX.Columns.Count Then
vertData = True
Else
vertData = False
End If
'Is the interpolation range verticle or horizontal
If interpX.Rows.Count > interpX.Columns.Count Then
vertInterp = True
Else
vertInterp = False
End If
'Is the X and Y area the same?
If dataRangeX.Rows.Count <> dataRangeY.Rows.Count And dataRangeX.Columns.Count <> dataRangeY.Columns.Count Then
SplineInterpolation = "#VALUE!"
Exit Function
End If
'Check for gradient variables
If IsMissing(gradient0) Then
'Natural spline interpolation
If vertData = True Then
numPoints = dataRangeX.Rows.Count - 1
ReDim Data(dataRangeX.Rows.Count)
ReDim h(dataRangeX.Rows.Count - 1)
ReDim alpha(dataRangeX.Rows.Count - 1)
ReDim l(dataRangeX.Rows.Count)
ReDim mu(dataRangeX.Rows.Count)
ReDim z(dataRangeX.Rows.Count)
ReDim b(numPoints)
ReDim c(numPoints)
ReDim d(numPoints)
For i = 0 To numPoints
Set Data(i) = New Point2D
Data(i).SetX (dataRangeX(i + 1, 1).Value)
Data(i).SetY (dataRangeY(i + 1, 1).Value)
Next
Else
numPoints = dataRangeX.Columns.Count - 1
ReDim Data(dataRangeX.Columns.Count)
ReDim h(dataRangeX.Columns.Count - 1)
ReDim alpha(dataRangeX.Columns.Count - 1)
ReDim l(dataRangeX.Columns.Count)
ReDim mu(dataRangeX.Columns.Count)
ReDim z(dataRangeX.Columns.Count)
ReDim b(numPoints)
ReDim c(numPoints)
ReDim d(numPoints)
For i = 0 To numPoints
Set Data(i) = New Point2D
Data(i).SetX (dataRangeX(1, i + 1).Value)
Data(i).SetY (dataRangeY(1, i + 1).Value)
Next
End If
'Does x only increase?
For i = 0 To numPoints - 1
If Data(i).GetX() >= Data(i + 1).GetX() Then
'If not then exit
SplineInterpolation = "#VALUE!"
Exit Function
End If
Next
For i = 0 To numPoints - 1
h(i) = Data(i + 1).GetX() - Data(i).GetX()
Next
For i = 1 To numPoints - 1
alpha(i) = (3 / h(i)) * (Data(i + 1).GetY() - Data(i).GetY()) - (3 / h(i - 1)) * (Data(i).GetY() - Data(i - 1).GetY())
Next
l(0) = 1
mu(0) = 0
z(0) = 0
For i = 1 To numPoints - 1
l(i) = 2 * (Data(i + 1).GetX() - Data(i - 1).GetX()) - h(i - 1) * mu(i - 1)
mu(i) = h(i) / l(i)
z(i) = (alpha(i) - h(i - 1) * z(i - 1)) / l(i)
Next
l(numPoints) = 1
z(numPoints) = 0
c(numPoints) = 0
For j = numPoints - 1 To 0 Step -1
c(j) = z(j) - mu(j) * c(j + 1)
b(j) = (Data(j + 1).GetY() - Data(j).GetY()) / h(j) - h(j) * (c(j + 1) + 2 * c(j)) / 3
d(j) = (c(j + 1) - c(j)) / (3 * h(j))
Next
Else
'Clamped spline interpolation
If vertData = True Then
numPoints = dataRangeX.Rows.Count - 1
ReDim Data(dataRangeX.Rows.Count)
ReDim h(dataRangeX.Rows.Count - 1)
ReDim alpha(dataRangeX.Rows.Count)
ReDim l(dataRangeX.Rows.Count)
ReDim mu(dataRangeX.Rows.Count)
ReDim z(dataRangeX.Rows.Count)
ReDim b(numPoints)
ReDim c(numPoints)
ReDim d(numPoints)
For i = 0 To numPoints
Set Data(i) = New Point2D
Data(i).SetX (dataRangeX(i + 1, 1).Value)
Data(i).SetY (dataRangeY(i + 1, 1).Value)
Next
Else
numPoints = dataRangeX.Columns.Count - 1
ReDim Data(dataRangeX.Columns.Count)
ReDim h(dataRangeX.Columns.Count - 1)
ReDim alpha(dataRangeX.Columns.Count)
ReDim l(dataRangeX.Columns.Count)
ReDim mu(dataRangeX.Columns.Count)
ReDim z(dataRangeX.Columns.Count)
ReDim b(numPoints)
ReDim c(numPoints)
ReDim d(numPoints)
For i = 0 To numPoints
Set Data(i) = New Point2D
Data(i).SetX (dataRangeX(1, i + 1).Value)
Data(i).SetY (dataRangeY(1, i + 1).Value)
Next
End If
'Does x only increase?
For i = 0 To numPoints - 1
If Data(i).GetX() >= Data(i + 1).GetX() Then
'If not then exit
SplineInterpolation = "#VALUE!"
Exit Function
End If
Next
fpo = gradient0
fpn = gradientN
For i = 0 To numPoints - 1
h(i) = Data(i + 1).GetX() - Data(i).GetX()
Next
alpha(0) = 3 * (Data(1).GetY() - Data(0).GetY()) / h(0) - 3 * fpo
alpha(numPoints) = 3 * fpn - 3 * (Data(numPoints).GetY() - Data(numPoints - 1).GetY()) / h(numPoints - 1)
For i = 1 To numPoints - 1
alpha(i) = (3 / h(i)) * (Data(i + 1).GetY() - Data(i).GetY()) - (3 / h(i - 1)) * (Data(i).GetY() - Data(i - 1).GetY())
Next
l(0) = 2 * h(0)
mu(0) = 0.5
z(0) = alpha(0) / l(0)
For i = 1 To numPoints - 1
l(i) = 2 * (Data(i + 1).GetX() - Data(i - 1).GetX()) - h(i - 1) * mu(i - 1)
mu(i) = h(i) / l(i)
z(i) = (alpha(i) - h(i - 1) * z(i - 1)) / l(i)
Next
l(numPoints) = h(numPoints - 1) * (2 - mu(numPoints - 1))
z(numPoints) = (alpha(numPoints) - h(numPoints - 1) * z(numPoints - 1)) / l(numPoints)
c(numPoints) = z(numPoints)
For j = numPoints - 1 To 0 Step -1
c(j) = z(j) - mu(j) * c(j + 1)
b(j) = (Data(j + 1).GetY() - Data(j).GetY()) / h(j) - h(j) * (c(j + 1) + 2 * c(j)) / 3
d(j) = (c(j + 1) - c(j)) / (3 * h(j))
Next
End If
'Create output data
If vertInterp = True Then
For i = 1 To interpX.Rows.Count
For j = 0 To numPoints - 1
If interpX(i, 1).Value >= Data(j).GetX() Then
x1 = Data(j).GetX()
y1 = Data(j).GetY()
jStop = j
'Exit For
End If
Next
If vertOutput = True Then
output(i, 1) = y1 + b(jStop) * (interpX(i, 1).Value - x1) + c(jStop) * (interpX(i, 1).Value - x1) ^ 2 + d(jStop) * (interpX(i, 1).Value - x1) ^ 3
Else
output(1, i) = y1 + b(jStop) * (interpX(i, 1).Value - x1) + c(jStop) * (interpX(i, 1).Value - x1) ^ 2 + d(jStop) * (interpX(i, 1).Value - x1) ^ 3
End If
Next
Else
For i = 1 To interpX.Columns.Count
For j = 0 To numPoints - 1
If interpX(1, i).Value >= Data(j).GetX() Then
x1 = Data(j).GetX()
y1 = Data(j).GetY()
jStop = j
'Exit For
End If
Next
If vertOutput = True Then
output(i, 1) = y1 + b(jStop) * (interpX(1, i).Value - x1) + c(jStop) * (interpX(1, i).Value - x1) ^ 2 + d(jStop) * (interpX(1, i).Value - x1) ^ 3
Else
output(1, i) = y1 + b(jStop) * (interpX(1, i).Value - x1) + c(jStop) * (interpX(1, i).Value - x1) ^ 2 + d(jStop) * (interpX(1, i).Value - x1) ^ 3
End If
Next
End If
SplineInterpolation = output
End Function
An Example
Using these points:
| X | Y |
| 1 | 0 |
| 2 | 3 |
| 3 | 2 |
| 4 | 7 |
| 5 | 5 |
| 6 | 3 |
and interpolating for these points:
| X | Interpolation |
| 1 | 0 |
| 1.3 | 1.35 |
| 1.7 | 2.68 |
| 2 | 3 |
| 2.3 | 2.59 |
| 2.7 | 1.84 |
| 3 | 2 |
| 3.3 | 3.29 |
| 3.7 | 5.71 |
| 4 | 7 |
| 4.3 | 7.135 |
| 4.7 | 6.07 |
| 5 | 5 |
| 5.3 | 4.17 |
| 5.7 | 3.43 |
gave the interpolation column.
Plotting this in Excel and cheating slightly by adding a smoothed line gives:
Excel interpolates slightly differently, although it looks like it uses splines. Adding more points to the interpolated points column would have made the line smoother without having to add the line to the chart.
If there are any errors, omissions, improvements that can be made, etc. then please comment.
Update
I forgot to post details of the custom class I made to hold details of the x,y data.
It's not particularly necessary to have done it this way but I was re-using a class I'd created for a clustering macro.
Anyway, here is the code:
Option Explicit
Private x As Double
Private y As Double
Private cluster As Integer
Public Sub SetX(temp As Double)
x = temp
End Sub
Public Sub SetY(temp As Double)
y = temp
End Sub
Public Function GetX() As Double
GetX = x
End Function
Public Function GetY() As Double
GetY = y
End Function
Public Sub SetC(temp As Integer)
cluster = temp
End Sub
Public Function GetC() As Integer
GetC = cluster
End Function
Apologies for any confusion.
Hi Sociolist,
ReplyDeletetried your code but there is a problem about the Dim Data() As Point2D
What kind of dim. is Point2D? Array, Integer?
Did you run this in Excel 2013 VBA?! Don't get it compiled...
Thanks
Sorry about that. I forgot to add the class I'd created to hold the x and y data. I've added it now and updated the introduction.
DeleteHi! This looks like exactly what I need for my purposes (interpolating from a calibration curve for scientific equipment). However, I'm having difficulty with the Point2D class. (I apologize - it's been over a decade since I've done anything with VBA, so my skills are, well, basic.) I've put the second set of code into an inserted class module, but it doesn't seem to accomplish defining the Point2D class. Is there some convention for naming the class module that I should know, or is there a line of code that needs to be added to define the class? Thanks for your work and your help. (BTW using Excel 2013 (the add-in to 2003 I'd used for this purpose no longer compiles)). Ian
ReplyDeleteI don't think so. It should be as simple as adding the class module and renaming to Point2D.
DeleteIf this doesn't work you can always replace the Point2D array data() with two double arrays datax() and datay(). I don't think it would make the code any more difficult. To be honest, it might make it a lot easier. Given the problems it seems to be causing I'll probably rewrite it without the custom class.