Too long to comment:
Here's another take, restricting to higher order splines. Suppose the polynomials comprised in the spline are given by $\{p_{0,1}(x),\cdots,p_{n-1,n}(x)\}$. Let $p_{i,i+1}(x) = a^{(0)}_{i,i+1} + a^{(1)}_{i,i+1}x + \cdots + a^{(n)}_{i,i+1}x^{n-1}$, $\forall ~i$, $n$ being odd. The variables here would, of course, be the coefficients of these polynomials (degree, we'll talk about it later). The following points are in order now:
(1) Interpolation constraints $\{x_i,f_i\}$ imply linear equality constraints on the coefficients. In particular those will be:
$$
a^{(0)}_{i,i+1} + a^{(1)}_{i,i+1}x_i + \cdots + a^{(n)}_{i,i+1}x_i^{n-1} = f_i, ~\forall ~i\\
a^{(0)}_{i,i+1} + a^{(1)}_{i,i+1}x_{i+1} + \cdots + a^{(n)}_{i,i+1}x_{i+1}^{n-1} = f_{i+1},~\forall ~i.
$$
$C^2$ requirement will imply another bunch of linear constraints. In particular those will be the following set of equations:
$$
a^{(1)}_{i-1,i} + 2a^{(2)}_{i-1,i}x_i + \cdots + (n-1)a^{(n)}_{i-1,i}x_i^{n-2} =
a^{(1)}_{i,i+1} + 2a^{(2)}_{i,i+1}x_i\cdots + (n-1)a^{(n)}_{i,i+1}x_{i}^{n-3},~\forall ~i,
$$
and
$$
2a^{(2)}_{i-1,i} + 6a^{(3)}_{i-1,i}x_i+ \cdots + (n-1)(n-2)a^{(n)}_{i-1,i}x_i^{n-3} =
2a^{(2)}_{i,i+1} + 6a^{(3)}_{i,i+1}+ \cdots + (n-1)(n-2)a^{(n)}_{i,i+1}x_{i}^{n-2},~\forall ~i.
$$
(2) Requirement of convexity over $[x_i,x_{i+1}]$, implies that the polynomial given by $p_{i,i+1}''(x)$ is positive over $[x_i,x_{i+1}]$. A uni-variate polynomial over $[x_i,x_{i+1}]$ is non-negative if and only if (see Victorial Powers and Bruce Reznick, Polynomials that Are Positive on an Interval for a nice explanation):
$$
p_{i,i+1}''(x) = g_{i,i+1}(x) + (x_{i+1}-x)(x-x_i)h_{i,i+1}(x),
$$
where $g_{i,i+1}$ and $h_{i,i+1}$ are SOS polynomials of degree at most $n$ and $n-2$, respectively. Now, $g_{i,i+1}(x)=z^\top G_{i,i+1}z$, where $G\succeq 0$ and $z=[1 ~x \cdots x^{(n-1)/2}]^\top$. Similarly, $h_{i,i+1}(x)=y^\top H_{i,i+1}y$, where $H_{i,i+1}\succeq 0$ and $y=[1 ~x \cdots x^{(n-1)/2-1}]^\top$. The matrices $G_{i,i+1}$ and $H_{i,i+1}$ are additional variables. Comparing coefficients on either side of the equation given by:
$$
2a^{(2)}_{i,i+1} + 6a^{(3)}_{i,i+1}+ \cdots + (n-1)(n-2)a^{(n)}_{i,i+1}x_{i}^{n-2} = g(x) + (x_{i+1}-x)(x-x_i)h(x),
$$
yields linear affine equations in coefficients of $p_{i,i+1}(x)$, $G$ and $H$. For brevity, let these equations be given as $\mathcal{L}_{i,i+1}\left(a^{(0)}_{i,i+1},\cdots,a^{(n)}_{i,i+1},G_{i,i+1},H_{i,i+1}\right)=0$. Note that you get a set of linear equations $\forall ~i$.
(3) An SDP solver (CVXPY, or alike) can now be used to find a feasible solution, i.e.:
$$
\min ~~1~~\mbox{subject to}\\
a^{(0)}_{i,i+1} + a^{(1)}_{i,i+1}x_i + \cdots + a^{(n)}_{i,i+1}x_i^{n-1} = f_i, ~\forall ~i\\
a^{(0)}_{i,i+1} + a^{(1)}_{i,i+1}x_{i+1} + \cdots + a^{(n)}_{i,i+1}x_{i+1}^{n-1} = f_{i+1},~\forall ~i\\
a^{(1)}_{i-1,i} + \cdots + (n-1)a^{(n)}_{i-1,i}x_i^{n-2} =
a^{(1)}_{i,i+1} + \cdots + (n-1)a^{(n)}_{i,i+1}x_{i}^{n-3},~\forall ~i\\
2a^{(2)}_{i-1,i} + \cdots + (n-1)(n-2)a^{(n)}_{i-1,i}x_i^{n-3} =
2a^{(2)}_{i,i+1} + \cdots + (n-1)(n-2)a^{(n)}_{i,i+1}x_{i}^{n-2},~\forall ~i\\
\mathcal{L}_{i,i+1}\left(a^{(0)}_{i,i+1},\cdots,a^{(n)}_{i,i+1},G_{i,i+1},H_{i,i+1}\right)=0\\
G_{i,i+1},H_{i,i+1} \succeq 0, \forall i.
$$
Lack of a feasible point, however, does not mean there does not exist a $C^2$ function, of a spline with degree higher than the one chosen or outside the realm of splines.
(4) The choice of the degree of the polynomial is something I don't have a clear answer for. I would assume that a degree greater than 8 should work -- for every polynomial piece, there are 2 interpolation constraints, 4 for $C^2$, and 2 LMI constraint.
(5) A disadvantage of this method is that its not practical to work with several thousand data points.
Hope it helps.
