if the functions to be "splined" are invertible on the intervals $[x_i,x_{i+1}]$ and, if $f_{i-1}(x)$ and $f_i(x)$ are two such functions that are invertible on $[x_{i-1},x_i]$, resp. on $[x_i,x_{i+1}]$, additionally $f_{i-1}^{-1}(x_i)=f_i^{-1}(x_i)$

then $$S(x)\ :=\ \mathcal{F}^{-1}\Big(\mathcal{S}\big(x_i,\mathcal{F}\left(y_i+\gamma(x_i)\,\right)\big)\Big)-\gamma(x) $$
yields an entire class of interpolating non-polynomial splines. I used the "strange" notation in order to stress the analogy to signal processing: if the original interpolation task is non-polynomial in the "time"-domain then "transforming" the $y_i$ via the inverse functions makes the interpolation polynomial in the "frequency" domain; interpreting interpolation as filtering completes the analogy to signal processing.

The key concept is Homomorphism and was introduced to signal processing in R.W. Schafer "Echo removal by discrete generalized linear filtering". Res. Lab. Electron. MIT, Tech. Rep. (1969)

The $\gamma(x)$ is a way to allow for mixed interpolation, e.g. exponentials plus polynomials; the case $\gamma(x)\equiv \mathrm{const}$ seem interesting in its own right when investigating the spline's limit behavior for e.g. $\gamma(x)\equiv \mathrm{const}\,\to\,\infty$, especially in the case of exponential splines.

Maybe its worth mentioning albeit trivial that this "homomorphic" splining directly carries over higher dimensions or parametric interpolation and one isn't limited to polynomials in the "frequency" domain; rational interpolants would be the "next bigger thing" to use.

**Addendum:**

the notation I used is aimed at emphasizing the analogy to transformation from the time domain to the frequency domain e.g. via a Fourier transform $\mathcal{F}$; the analogy to interpolating $(x_i,y_i)$ with exponential splines $S(x)\in C^{k-1}$ of the form $e^{\sum_{j=0}^k a_{ij}x^j}, S(x_i)=y_i$ is that the tranformation $\mathcal{F}$ corresponds to applying the analogy of the Fourier transformation $\mathcal{F}$, i.e. the inverse $\ln(\cdot)$ of $e^{(\cdot)}$, to the ordinates $y_i$ of the data to be interpolated and then calculate the interpolating polynomial spline with these transformed ordinates; that spline is denoted by $\mathcal{S}(x)$ because it is calculated in the analogy of the frequency domain.

the final step is to go back to the analogy of the time domain by applying the analogy $\mathcal{F}^{-1}$, i.e. $e^{(\cdot)}$ of the inverse Fourier transformation to $\mathcal{S}(x)$ to obtain the interpolating spline $S(x)$ in the analogy of the time domain.

The $\gamma(x)$ allow for modeling situations where the data is "almost" polynomial, i.e. if for some order $h$ the magnitudes of the divided differences of order $h$ become small, then it may, depending on the model of the origin of the data, improve the quality of the interpolating spline if the values of a "lower polynomial hull" at the abscissas $x_i$ are subtracted from the ordinates prior to the transition to the analogy of the frequency domain.

Another usecase for introducing $\gamma(x)$ is if one is only interested in the shape of the interpolating curve but the ordinates of the data are outside the value range of the interpolation curve, e.g. negative $y_i$ in the case of exponential splines; then adding $\gamma(x)$ can fix these issues.

An example may make things clearer:

assuming we have sampled an empirical distribution of which we assume that is the sum of Gaussian's that we would like to recover.
In that case an Ansatz could be to take the logarithm of the ordinates, calculate a quadratic spline and take the abscissas of apices of parabolas with negative leading coefficient as an initial guess for the central values of the Gaussians.

2more comments