Cubic spline interpolation - how to calculate second derivative

I ask this qeustion on stackexchange sites: stackoverflow, codereview, and signal processing and no one can help and they send me here :)

So I implement cubic spilne interpolation in Java base on this document: http://www.geos.ed.ac.uk/~yliu23/docs/lect_spline.pdf

At first author shows how to calculate linear spline interpolation and I did this on my data, and receive this result:

It should be similar to this:

The overall shape is good but to receive better results I should use cubic spilne intepolation with is extend of linear interpolation but here problems starts. I don't understand how they calculate second derivative.

Could you explain that to me?

I also implement cubic spilne interpolation based on algorithm from wiki: http://en.wikipedia.org/w/index.php?title=Spline_%28mathematics%29&oldid=288288033#Algorithm_for_computing_natural_cubic_splines

x(0) is points.get(0), y(0) is values[points.get(0)], α is alfa and μ is mi Rest is the same as in pseudocode.

    public void createSpline(double[] values, ArrayList<Integer> points){
a = new double[points.size()+1];

for (int i=0; i <points.size();i++)
{
a[i] = values[points.get(i)];

}

b = new double[points.size()];
d = new double[points.size()];
h = new double[points.size()];

for (int i=0; i<points.size()-1; i++){
h[i] = points.get(i+1) - points.get(i);
}

alfa = new double[points.size()];

for (int i=1; i <points.size()-1; i++){
alfa[i] = (double)3 / h[i] * (a[i+1] - a[i]) - (double)3 / h[i-1] *(a[i+1] - a[i]);
}

c = new double[points.size()+1];
l = new double[points.size()+1];
mi = new double[points.size()+1];
z = new double[points.size()+1];

l[0] = 1; mi[0] = z[0] = 0;

for (int i =1; i<points.size()-1;i++){
l[i] = 2 * (points.get(i+1) - points.get(i-1)) - h[i-1]*mi[i-1];
mi[i] = h[i]/l[i];
z[i] = (alfa[i] - h[i-1]*z[i-1])/l[i];
}

l[points.size()] = 1;
z[points.size()] = c[points.size()] = 0;

for (int j=points.size()-1; j >0; j--)
{
c[j] = z[j] - mi[j]*c[j-1];
b[j] = (a[j+1]-a[j]) - (h[j] * (c[j+1] + 2*c[j])/(double)3) ;
d[j] = (c[j+1]-c[j])/((double)3*h[j]);
}

for (int i=0; i<points.size()-1;i++){
for (int j = points.get(i); j<points.get(i+1);j++){
//                fk[j] = values[points.get(i)];
functionResult[j] = a[i] + b[i] * (j - points.get(i))
+ c[i] * Math.pow((j - points.get(i)),2)
+ d[i] * Math.pow((j - points.get(i)),3);
}
}

}


but result look like (That's the other part of the signal) :

when is should be similar to this:

• And as I already told you on your other post, your wikipedia reference is outdated and of dubious quality, as even stated on that page. If you want to compute an approximative spline, which would be sensible for this data, then you have to look for that concept. – LutzL Dec 18 '13 at 13:56
• See (mathworks.de/de/help/curvefit/…) for a short introduction. – LutzL Dec 18 '13 at 14:08

According to the paper, which is a much much better reference than the outdated version of the wikipedia article, you got to solve

$$\frac{x_{k}-x_{k-1}}6 a_{k-1}+\frac{x_{k+1}-x_{k-1}}3 a_k + \frac{x_{k+1}-x_k}6 a_{k+1}=\frac{y_{k}-y_{k-1}}{x_{k}-x_{k-1}}-\frac{y_{k+1}-y_k}{x_{k+1}-x_k}$$

for $k=1,...,N-1$ and set $a_0=0=a_N$. To this end, you can implement the Jacobi or Gauss-Seidel method for linear systems.

$$a_k^{new} = \frac3{x_{k+1}-x_{k-1}}\left(\frac{y_{k}-y_{k-1}}{x_{k}-x_{k-1}}-\frac{y_{k+1}-y_k}{x_{k+1}-x_k}-\frac{x_{k}-x_{k-1}}6 a_{k-1}-\frac{x_{k+1}-x_k}6 a_{k+1}\right)$$

where you either compute all new values before replacing the current vector of $a$, or you can set $a_k=a_k^{new}$ directly after its computation.

Correction: There was a sign error on the right side, it must be the difference of the slopes; the paper has it correctly.

Remark: The system is visibly diagonal dominant, convergence should be fast.

Added: How to get a spline approximation: First variant: Use B-splines as basis functions. Second variant: Use the parametrization

$$a+bx+cx^2+\sum_{k=0}^M d_k\,(x-x_k)_+^3,\quad \text{where }(x-x_k)_+=\max(0,x-x_k)$$

with a sparse distribution of $x_0,...,x_M$ relative to the amount of sample points, and then do a standard least squares fit. If $x_M$ is outside the interval of sample points, then leave out the last term and set $d_M=-\sum_{k=0}^{M-1}d_k$ to get a proper spline with quadratic ends.

There is a simple convolutional method for generating a natural cubic spline interpolation. See:

http://nbviewer.jupyter.org/github/mtimmerm/IPythonNotebooks/blob/master/NaturalCubicSplines.ipynb