Compute Fourier coefficients of spline fit to data

Question

Suppose you have data $$\{(x_i, -1^{i+1})\}_{1\dots N}, \quad x_1=0<x_i<x_{i+1}<x_N=2\pi \ \forall i \in\{2,\dots N-1\} $$ In other words, we have a sequence of $y=\pm1$ values at distinct random locations in $[0,2\pi]$, with the points $P_1=(0,1)$ and $P_N=(2\pi,-1^{N+1})$ always included in the sequence. This data set is interpolated with splines, specifically MATLAB pchip splines. The interpolation spline so obtained is periodic, thus it makes sense to compute its Fourier coefficients. I have been asked to perform such a computation, and I was thinking to use MATLAB fft function to do this, after sampling the spline interpolant at a "large enough" (?) number of points. Questions:

• the pchip interpolant is only $C^1$ (unlike cubicspline which is $C^2$), thus the Fourier series will converge more slowly. Should I suggest that cubicspline be used for the interpolation step, instead of pchip?
• is using fft really necessary? After all, interpolating splines have a relatively simple analytical expression. I was wondering if there could be an explicit formula for the computation of the Fourier coefficients... Maybe it could be possible to arrive at some linear system? Hopefully, solving it by A\b could be faster and/or more accurate than fft, if the dimension of A is small enough.

Answer 1

Ok maybe I understand a bit now. I will make a try. You don't have to use fft. Consider the following:

$$f(t) = \sum_{\forall k} \left(c_k\cos(kt)+s_k\sin(kt)\right)$$

Now if we consider the requirements on the individual time points $t_l$:

$$f(t_l) = \left\{\pm 1\right\}_{l} \Leftrightarrow \sum_{\forall k} \left(\underset{\text{Constant for each }t_l,k}{c_k\cdot\underbrace{\cos(kt_l)}} + \underset{\text{Constant for each }t_l,k}{s_k\cdot\underbrace{\sin(kt_l)}}\right) = {\{\pm 1\}}_l$$

Calculate the amplitudes of each sine / cosine at each position you will be getting a linear equation system with the coefficients times the amplitudes of the wave at each position summed - one linear equation per spatial point. This should be exactly solvable when you add as many overtones so you get as many equations as $\pm 1$ positions. So you can write it as $\bf Ax=b$ where 1 row in $\bf A$ for each $\pm 1$ impulse and two columns for each new overtone (one sin and one cos). You can probably calculate for how few number of overtones you would require to solve it, but as you mention you may want some kind of a regularizer.

Here is a plot without any regularization:

So we see we may need some regularization of some kind...

However, differentiation is easy to perform in fourier domain. Maybe you can set first derivative = 0 at the $\pm 1$ points. That would give a new set of equations which you probably can derive quite easily from above expression. It is probably easier to encourage derivative to be 0 ( to create local maxima and minima ) than to try and do a fitting to FFT of piecewise polynomials.

---

If this is still not enough, then maybe put blanket second order derivatives should be close to 0 everywhere except at the maxima / minima points of interest - or maybe only in the middle points of the $\pm 1$ points. Same as previous idea, since differentiation is an eigenoperator on our basis functions, this will remain linear, it will just grow the equation system a bit.

---

Fourier transforms of discontinuous or discrete functions will be jumpy due to the Gibbs phenomenon. Yet another improvement would be to switch kernel at some point of the processing. Fejér kernels are for example much smoother. You can see some other weighted kernels I developed in this question.

---

What would be nasty about expressing piecewise polynomials in the Fourier domain is that the multiplication in the monomials $t^n$ will turn into convolution in the Fourier domain even before you start windowing (multiplying with sinc) or try to build linear combinations of them. It will all become very convoluted and non-linear.