Loading web-font TeX/Math/Italic

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