From 9bd5546470b050adf68c5109e98fe95220967d5b Mon Sep 17 00:00:00 2001
From: Ahmet Inan <inan@aicodix.de>
Date: Wed, 28 Mar 2018 16:29:46 +0200
Subject: [PATCH] added Algorithm for computing Natural Cubic Splines

Modified algorithm from Wikipedia to work with integer x_i and x_i >= 0:
https://en.wikipedia.org/wiki/Spline_(mathematics)#Algorithm_for_computing_natural_cubic_splines
---
 README.md |  5 +++++
 spline.hh | 57 +++++++++++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 62 insertions(+)
 create mode 100644 spline.hh

diff --git a/README.md b/README.md
index d6ef8ec..38342df 100644
--- a/README.md
+++ b/README.md
@@ -33,3 +33,8 @@ Interface for reading and writing [PCM](https://en.wikipedia.org/wiki/Pulse-code
 
 Read and write [WAV](https://en.wikipedia.org/wiki/WAV) files
 
+### [spline.hh](spline.hh)
+
+Algorithm for computing uniform and [natural cubic splines](https://en.wikipedia.org/wiki/Spline_(mathematics)#Algorithm_for_computing_natural_cubic_splines)
+Very useful for data interpolation.
+
diff --git a/spline.hh b/spline.hh
new file mode 100644
index 0000000..1093d68
--- /dev/null
+++ b/spline.hh
@@ -0,0 +1,57 @@
+/*
+Some spline algorithms
+
+Copyright 2018 Ahmet Inan <inan@aicodix.de>
+*/
+
+#ifndef SPLINE_HH
+#define SPLINE_HH
+
+namespace DSP {
+
+template <int KNOTS, typename OTYPE, typename ITYPE>
+class UniformNaturalCubicSpline
+{
+	OTYPE A[KNOTS], B[KNOTS], C[KNOTS], D[KNOTS];
+public:
+	UniformNaturalCubicSpline(OTYPE *y, int STRIDE = 1)
+	{
+		for (int i = 0; i < KNOTS; ++i)
+			A[i] = y[i * STRIDE];
+		ITYPE u[KNOTS], l[KNOTS];
+		u[0] = 0;
+		l[0] = l[KNOTS-1] = 1;
+		OTYPE z[KNOTS];
+		z[0] = z[KNOTS-1] = 0;
+		for (int i = 1; i < KNOTS - 1; ++i) {
+			l[i] = ITYPE(4) - u[i-1];
+			u[i] = ITYPE(1) / l[i];
+			z[i] = (ITYPE(3) * (A[i+1] - ITYPE(2) * A[i] + A[i-1]) - z[i-1]) / l[i];
+		}
+		C[KNOTS-1] = 0;
+		for (int i = KNOTS - 2; i >= 0; --i) {
+			C[i] = z[i] - u[i] * C[i+1];
+			B[i] = A[i+1] - A[i] - (C[i+1] + ITYPE(2) * C[i]) / ITYPE(3);
+			D[i] = (C[i+1] - C[i]) / ITYPE(3);
+		}
+	}
+	OTYPE operator () (ITYPE x)
+	{
+		int k = x;
+		ITYPE t = x - ITYPE(k);
+		if (k < 0) {
+			t = x;
+			k = 0;
+		}
+		if (k >= KNOTS - 1) {
+			t = x - ITYPE(KNOTS-2);
+			k = KNOTS-2;
+		}
+		return A[k] + t * (B[k] + t * (C[k] + t * D[k]));
+	}
+};
+
+}
+
+#endif
+