# Compute Fourier coefficients of spline fit to data

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.
• Most of what you say, seems reasonable to me. It is a good question whether "oversampling" is the way to go here, and will depend on the application (target accuracy vs work; remember, however, that fft is fast). You could compute the Fourier coefficient analytically, of course: after all, the input is a piecewise polynomial. – user66081 Jul 22 '16 at 19:41
• @user66081, yes, I'd like to do that but it's been quite some since I last computed Fourier coefficients analytically. The spline interpolant is a linear combination of spline basis functions, which locally are third degree polynomials, but globally on $[0,2\pi]$ are "something else"...could you write an answer showing how to compute the Fourier coefficients analytically? – DeltaIV Jul 22 '16 at 21:58
• @user66081,I'm not sure how it is relevant. Sinc interpolation works for bandlimited signals. My signal is a power signal, not an energy signal, and it's periodic: that's why I'm talking about Fourier series and not Fourier transform. Am I missing something? – DeltaIV Jul 24 '16 at 11:26
• Do you just want to inverse-transform or what is the objective? If it just is a Fourier series it is just a linear combination of sines and cosines so you can just generate them and weigh them together with the $\pm 1$. – mathreadler Jul 29 '16 at 7:47

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.

• Nice! Expecially the point on adding 0 first derivative equations at the $\pm1$ points. I don't understand the part about kernels though - remember that my objective is to compute spectra. How do I use your kernels to do that, i.e., how do I modify $bf Ax=b$ to use kernels? – DeltaIV Jul 31 '16 at 10:11
• The vectorization is essentially with respect to a vector with Fourier coefficients. You can build matrices which transform or move your coefficients as you want. For example a DFT matrix or a set of FFT factor matrices or their inverses or some third thing. If you switch to solving the normal equations $A^TAx = A^Tb$ you can add other constraints as you see fit but it may get complicated if you are not used thinking like that. But if you are OK to think like that you could probably add polynomial fitting in another domain into the same equation system. – mathreadler Jul 31 '16 at 10:20
• The kernels are how you choose to transform the spectra back into the spatial or temporal domain. Like a reweighting of the Fourier coefficients after you have found the spectra to get nicer behaviour in the ordinary domain. To avoid effects like ringing and gibbs phenomenon - but at a cost of worse square integral fit. – mathreadler Jul 31 '16 at 10:41
• Thank you! One last question: since the linear combination of sines and cosines is always periodic, should I discard the last equation, i.e., the one which corresponds to passage through point $(2\pi,1)$, since I already have an equation which corresponds to passage through point $(0,1)$? – DeltaIV Jul 31 '16 at 13:03
• Yes as long as you use fourier series of a compatible base frequency you should be OK doing that. But when you start to introduce regularizers so you may not be sure to get exact fit at all points any longe. Then the extra point will make a good fit in the seams have double importance as compared to the other points. – mathreadler Jul 31 '16 at 13:08