Here is also an interpolator ported from C++ to Java with the original version in the book “Numerical Recipies”:
public abstract class BaseInterpolator {
int n;
int mm;
int jsav;
int cor;
int dj;
double[] xx;
public BaseInterpolator(double[] x, int m) {
n = x.length;
mm = m;
jsav = 0;
cor = 0;
xx = x;
dj = Math.min(1, (int) Math.pow((double) n, 0.25));
}
public double interpolate(double x) {
int jlo = (cor != 0) ? hunt(x) : locate(x);
return rawInterpolate(jlo, x);
}
public void interpolate(double x, double[] input, double[] result) {
int jlo = (cor != 0) ? hunt(x) : locate(x);
rawInterpolate(jlo, x, input, result);
}
public double interpolate(double x, double[] input) {
int jlo = (cor != 0) ? hunt(x) : locate(x);
return rawInterpolate(jlo, x, input);
}
public int locate(final double x) {
int ju, jm, jl;
if (n < 2 || mm < 2 || mm > n) {
System.out.println(this + " Locate size error");
return 0;
}
boolean ascnd = (xx[n - 1] >= xx[0]);
jl = 0;
ju = n - 1;
while (ju - jl > 1) {
jm = (ju + jl) >> 1;
if (x >= xx[jm] == ascnd) {
jl = jm;
} else {
ju = jm;
}
}
cor = Math.abs(jl - jsav) > dj ? 0 : 1;
jsav = jl;
return Math.max(0, Math.min(n - mm, jl - ((mm - 2) >> 1)));
}
public int hunt(final double x) {
int jl = jsav;
int jm, ju;
int inc = 1;
if (n < 2 || mm < 2 || mm > n) {
System.out.println(this + " Hunt size error");
return 0;
}
boolean ascnd = (xx[n - 1] > xx[0]);
if (jl < 0 || jl > n - 1) {
jl = 0;
ju = n - 1;
} else {
if (x >= xx[jl] == ascnd) {
for (;;) {
ju = jl + inc;
if (ju >= n - 1) {
ju = n - 1;
break;
} else if (x < xx[ju] == ascnd) {
break;
} else {
jl = ju;
inc += inc;
}
}
} else {
ju = jl;
for (;;) {
jl = jl - inc;
if (jl <= 0) {
jl = 0;
break;
} else if (x >= xx[jl] == ascnd) {
break;
} else {
ju = jl;
inc += inc;
}
}
}
}
while (ju - jl > 1) {
jm = (ju + jl) >> 1;
if (x >= xx[jm] == ascnd) {
jl = jm;
} else {
ju = jm;
}
}
cor = Math.abs(jl - jsav) > dj ? 0 : 1;
jsav = jl;
return Math.max(0, Math.min(n - mm, jl - ((mm - 2) >> 1)));
}
public double rawInterpolate(int jlo, double x) {
return 0.0;
}
public double rawInterpolate(int jlo, double x, double[] input) {
return 0.0;
}
public void rawInterpolate(int jlo, double x, double[] input,
double[] result) {
}
}
public abstract class DoubleBaseInterpolator extends BaseInterpolator {
double[] yy;
public DoubleBaseInterpolator(double[] x, double[] y, int m) {
super(x, m);
this.yy = y;
}
}
public class CubicSplineInterpolator extends DoubleBaseInterpolator {
double[] y2;
public CubicSplineInterpolator(double[] xValues, double[] yValues,
double yp1, double ypn) {
super(xValues, yValues, 2);
y2 = new double[xValues.length];
sety2(xValues, yValues, yp1, ypn);
}
public CubicSplineInterpolator(double[] xValues, double[] yValues) {
super(xValues, yValues, 2);
y2 = new double[xValues.length];
sety2(xValues, yValues, 1.e99, 1.e99);
}
private void sety2(double[] xv, double[] yv, double yp1, double ypn) {
int i, k;
double p, qn, sig, un;
int n = y2.length;
double[] u = new double[n - 1];
if (yp1 > 0.99e99) {
y2[0] = u[0] = 0.0;
} else {
y2[0] = -0.5;
u[0] = (3.0 / (xv[1] - xv[0]))
* ((yv[1] - yv[0]) / (xv[1] - xv[0]) - yp1);
}
for (i = 1; i < n - 1; i++) {
sig = (xv[i] - xv[i - 1]) / (xv[i + 1] - xv[i - 1]);
p = sig * y2[i - 1] + 2.0;
y2[i] = (sig - 1.0) / p;
u[i] = (yv[i + 1] - yv[i]) / (xv[i + 1] - xv[i])
- (yv[i] - yv[i - 1]) / (xv[i] - xv[i - 1]);
u[i] = (6.0 * u[i] / (xv[i + 1] - xv[i - 1]) - sig * u[i - 1]) / p;
}
if (ypn > 0.99e99) {
qn = un = 0.0;
} else {
qn = 0.5;
un = (3.0 / (xv[n - 1] - xv[n - 2]))
* (ypn - (yv[n - 1] - yv[n - 2]) / (xv[n - 1] - xv[n - 2]));
}
y2[n - 1] = (un - qn * u[n - 2]) / (qn * y2[n - 2] + 1.0);
for (k = n-2; k >= 0; k--) {
y2[k] = y2[k] * y2[k + 1] + u[k];
}
}
@Override
public double rawInterpolate(int jl, double x) {
int klo = jl, khi = jl + 1;
double y, h, b, a;
h = xx[khi] - xx[klo];
if (h == 0.0) {
System.out.println(this + " bad input");
return 0.0;
}
a = (xx[khi] - x) / h;
b = (x - xx[klo]) / h;
y = a * yy[klo] + b * yy[khi]
+ ((a * a * a - a) * y2[klo] + (b * b * b - b) * y2[khi])
* (h * h) / 6.0;
return y;
}
public static void main(String[] args) {
CubicSplineInterpolator li = new CubicSplineInterpolator(new double[] {
0.0, 1.2, 3.0, 5.0}, new double[] { 5.0, 2.0, -11.0, 7 });
for (int i = 0; i < 300; i++) {
double x = i * 0.01;
System.out.println("input: " + x + " output: " + li.interpolate(x));
}
}
}
The problem with this is that it is slower than the specialized version in my previoius post but it allows a non-even distribution of the “x-values”. The other is just a uniform, but very quick, variant.
