|
|
Spline interpolation for smooth display
In test and measurement
applications it is sometimes desirable to compute a smoothed curve from the raw measurement data. Also, we sometimes require information in between curve points. This can be achieved by applying an
interpolation algorithm to the raw curve data.
One of the most popular
methods is the 3rd order spline interpolation algorithm. This algorithm constructs segments of 3rd order polynomials which are continuous up to the second derivative and is rather well-behaved.
Please read the literature for more information on the theory behind this and other interpolation methods.
Functions
SplineCompute2nDerivative and Splint which you find below were ported from chapter 3 of "Numerical Recipes in C / Cambridge University Press".
Computing an interpolated version of a curve using these functions requires the following steps
Obtain nRaw points of the original data and place the curve points in arrays xRaw and yRaw.
From the original data
compute the array y2Raw of second derivatives of the original curve using function SplineCompute2nDerivative.
Set up a dense grid of nGrid x-values xGrid for which the interpolated function shall be evaluated.
For each abscissa xGrid evaluate the interpolated function value yGrid using function Splint.
Pass arrays xGrid and yGrid to DynaPlot.
In Visual Basic this could be coded like this
|
Const nRaw = 20
Const nGrid = 100
Dim xRaw(nRaw - 1) As Double
Dim yRaw(nRaw - 1) As Double
Dim y2Raw(nRaw - 1) As Double
Dim xGrid(nGrid - 1) As Double
Dim yGrid(nGrid - 1) As Double
Private Sub Form_Load()
Dim i As Long
Dim v As Double
v = nRaw / nGrid
Dim yp1 As Double ' the first derivative at the beginning of the curve
Dim ypn As Double ' the first derivative at the end of the curve
yp1 = 1E+30 ' compute natural spline for this boundary
ypn = 1E+30 ' compute natural spline for this boundary
' simulate some measured data
For i = 0 To nRaw - 1
xRaw(i) = i
yRaw(i) = 0.01 * i * i + Sin(i / 2) * Sin(i / 2)
Next
SplineCompute2nDerivative xRaw, yRaw, nRaw, yp1, ypn, y2Raw
' Compute the interpolated function on a dense grid
For i = 0 To nGrid - 1
xGrid(i) = i * v
yGrid(i) = Splint(xRaw, yRaw, y2Raw, nRaw, xGrid(i))
Next
' plot the functions
DynaPlot1.DataCurves.Add "Raw", xRaw, yRaw,-1, True
DynaPlot1.DataCurves.Add "Interpolated", xGrid, yGrid,-1, True
DynaPlot1.Legend.Show.State = True
End Sub
|
Functions SplineCompute2nDerivative is called once to obtain the second derivative of the original
function given by arrays x and y. Argument n is the length of these arrays and also the length of result
array y2. Arguments yp1 and ypn are the first derivatives of the interpolating function at the first and
last curve points respectively. A value >= 1e30 signals that a natural spline shall be used.
|
Public Sub
SplineCompute2nDerivative(x As Variant, y As Variant, n As Long, yp1 As
Double, ypn As Double, y2 As Variant)
Dim i As Long, k As Long
Dim p As Double, qn As Double, sig As Double, un As Double
ReDim u(n - 1) As Double
If yp1 > 9.9E+29 Then
y2(0) = 0#
u(0) = 0#
Else
y2(0) = -0.5
u(0) = (3# / (x(1) - x(0))) * ((y(1) - y(0)) / (x(1) - x(0)) - yp1)
End If
For i = 1 To n - 2
sig = (x(i) - x(i - 1)) / (x(i + 1) - x(i - 1))
p = sig * y2(i - 1) + 2#
y2(i) = (sig - 1#) / p
u(i) = (y(i + 1) - y(i)) / (x(i + 1) - x(i)) - (y(i) - y(i - 1)) /
(x(i) - x(i - 1))
u(i) = (6# * u(i) / (x(i + 1) - x(i - 1)) - sig * u(i - 1)) / p
Next
If ypn > 9.9E+29 Then
qn = 0
un = 0#
Else
qn = 0.5
un = (3# / (x(n) - x(n - 1))) * (ypn - (y(n) - y(n - 1)) / (x(n) - x(n - 1)))
End If
y2(n - 1) = (un - qn * u(n - 1)) / (qn * y2(n - 1) + 1#)
For k = n - 2 To 0 Step -1
y2(k) = y2(k) * y2(k + 1) + u(k)
Next
End Sub
|
Function Splint is called for every abscissa x of the dense grid. Arguments to the function are xa and ya,
the arrays of point coordinates of the original curve and n, the size of these arrays. The function returns the interpolated ordinate value at abscissa x.
|
Public Function
Splint(xa As Variant, ya As Variant, y2a As Variant, n As Long, ByVal x As Double) As Double
Dim klo As Long, khi As Long, k As Long
Dim h As Double, b As Double, a As Double
klo = 0
khi = n - 1
While khi - klo > 1
k = (khi + klo) / 2
If xa(k) > x Then khi = k Else: klo = k
Wend
h = xa(khi) - xa(klo)
If h = 0# Then MsgBox ("Bad xa input to routine Splint"): Exit Function
a = (xa(khi) - x) / h
b = (x - xa(klo)) / h
Splint = a * ya(klo) + b * ya(khi) + ((a * a * a - a) * y2a(klo) + (b * b *
b - b) * y2a(khi)) * (h * h) / 6#
End Function
|
The above code can also be downloaded:
Download VB demo project
|
