Thursday, 11 October 2012

Cubic splines: Interpolation and smoothing

Cubic splines are widely used technique for piecewise smooth interpolation. Between any two points of a discretely given function, the function is approximated by a third order polynomial. However, splines can also be used for smoothing or fitting of measurements. The theory behind the splines is nicely explained by Pollock (1999). I implemented that algorithm in IDL procedure splinecoeff.pro.

To try it yourself, here is a small example. Download the text file with the values of the sunspot area data, daily_area.txt from solarscience.msfc.nasa.gov, and run the following:

; Load the data
n = FILE_LINES('daily_area.txt')-1
OPENR, 1, 'daily_area.txt'
a = ''
READF, 1, a
area = FLTARR(6, n)
READF, 1, area
CLOSE, 1

; Extract data between 1977-01-01 and 2011-12-31
initial = WHERE(area[0, *] EQ 1977 AND area[1, *] EQ  1 AND area[2, *] EQ  1)
final   = WHERE(area[0, *] EQ 2011 AND area[1, *] EQ 12 AND area[2, *] EQ 31)
time = JULDAY(area[1,*], area[2,*], area[0,*])
x = time[initial:final]
y = REFORM(area[3, initial:final])
nx = N_ELEMENTS(x)

; Compute coefficients of the cubic splines
coeff = SPLINECOEFF(x, y, lambda = 1.d3)

; Plot the original data and the spline function for given lambda
y1 = FLTARR(nx-1)
x1 = x[0:nx-2]
FOR i = 0, N_ELEMENTS(y)-2 DO y1[i] = coeff.d[i] + $
                                      coeff.c[i] * (x[i+1]-x[i]) + $
                                      coeff.b[i] * (x[i+1]-x[i])^2 + $
                                      coeff.a[i] * (x[i+1]-x[i])^3
PLOT, x, y, psym = 3
OPLOT, x1, y1
The only difference between the interpolation, smoothing and fitting of the measurement is in the free parameter $\lambda$. If $\lambda$ is very small, than the result is the cubic spline interpolation. If $\lambda$ is very large, the measurement is smoothed by cubic splines. In the extreme case, the resulting splines become a linear function. The following animation show how the result change with $\lambda$:



Note: See paper of Djurovic and Paquet (1999) for an example where the cubic splines are used for filtering and smoothing of data.

No comments:

Post a Comment