trigonometric interpolation with non-equidistant sampling

I want to interpolate a periodic function on a non-equidistant grid and have implemented it the using the Lagrangian formula described in wikipedia. For an odd number of data points, this takes the form

$$p(x) = \sum_{k=0}^{2K} y_k\,t_k(x)$$

with

$$t_k(x) = \prod_{m=0,m\ne k}^{2K} \frac{\sin\frac12(x-x_m)}{\sin\frac12(x_k-x_m)}$$

which bears strong resemblance to the usual Lagrangian interpolation formula, but with sines and a the factors 1/2 appearing. This formula makes perfect sense and in fact, I used the equivalent of the barycentric form described in the wikipedia page of Lagrange polynomials (I cannot post two links):

$$p(x) = \frac{\sum_{k=0}^{2K} \frac{w_k}{\sin{\frac{1}{2}\left(x-x_k\right)}}y_k}{\sum_{k=0}^{2K} \frac{w_k}{\sin{\frac{1}{2}\left(x-x_k\right)}}}$$

with

$$w_k = \frac{1}{l'\left(x_k\right)}$$

and

$$l\left(x\right) = \prod_{k=0}^{2K} \sin{\frac{1}{2}\left(x-x_k\right)}$$ .

This works magnificently, as can be seen in the following figure, where I use 15 (red) points of

$$f\left(x\right) = \sin{\left(x\right)} + \frac{1}{4}\cos{\left(2x\right)} + 2 \sin{\left(6x\right)} + 3 \cos{\left(7 x\right)}$$ ,

defined on a non-equidistant grid (in this case with more points on the edges of the region $2\ldots \pi$). As Nyquist states, this is indeed enough information to interpolate $f\left(x\right)$, which I illustrate on 100 (blue) points, also non-equidistant (but in this case with more points in the middle of the region).

However, for an even number of points, there is an additional degree of freedom. In the Wikipedia entry for trigonometric interpolation, it is stated that this is usually employed to make sure the highest modal content only has a cosine, and no sine.

This, they state, is implemented by using additional factors

$$\alpha_k=\sum_{m=0,m\ne k}^{2K-1} x_m$$

that modify the formulas to

$$t_k(x) = \frac{\sin\frac12(x-\alpha_k)}{\sin\frac12(x_k-\alpha_k)}\prod_{m=0,m\ne k}^{2K-1} \frac{\sin\frac12(x-x_m)}{\sin\frac12(x_k-x_m)}$$

This, however, I have tried, and I have found it does not work. Not only does the interpolated function not go through the correct points, but also it blows up to high absolute values.

Therefore, I have investigated it in more detail, and have found in fact that this definition of $\alpha_k$ contains singularities.

Take for example the case with 4 points

$$x_k = \frac{k \pi}{4} \ , k = 0\ldots 3$$

applied to a simple sine. Note that according to Nyquist, 3 points would be enough, so there is indeed a redundancy of 1 point.

Let's again look at the interval $x= 0\ldots 2 \pi$. As the only non-vanishing components of $y_k$ are 1 and 3 (with $y_1 = 1$ and $y_3=-1$), the necessary factors $t_k$ are given by

$$t_1 = \sqrt{2} \sin{\frac{x}{2}} \cos{\frac{x}{2}} \left( \cos{\frac{x}{2}} + \sin{\frac{x}{2}} \right ) \frac{\sin{\frac{1}{2}\left(x-\alpha_1\right)}}{\sin{\frac{1}{2}\left(\frac{\pi}{2} - \alpha_1\right)}} \\ t_3 = \sqrt{2} \sin{\frac{x}{2}} \cos{\frac{x}{2}} \left( \cos{\frac{x}{2}} - \sin{\frac{x}{2}} \right ) \frac{\sin{\frac{1}{2}\left(x-\alpha_3\right)}}{\sin{\frac{1}{2}\left(\frac{3\pi}{2} - \alpha_3\right)}}$$

where it can be seen that there are extra degrees of freedom $\alpha_k$.

However, now using their definition:

$$\alpha_1 = \frac{5}{2} \pi \\ \alpha_3 = \frac{3}{2} \pi$$ ,

this clearly leads to singularities in the denominators of $t_1$ and $t_3$. Does anybody know why this happens? There must either be something wrong in the definition of $\alpha_k$ or in my understanding of it.

I have tried to find other resources but the only thing I could find limited the discussion to an odd number of points (which indeed works fine.)

I have also tried to follow the arguments in the wikipedia page and I am not able to capture why they reach that expression for $\alpha$, repeated here:

$$\alpha_k=\sum_{m=0,m\ne k}^{2K-1} x_m$$

By the way, as a side note, it can be seen that throwing away the fourth point $\frac{3}{2}\pi$, the remaining three points indeed lead to the correct interpolation:

$$p\left(x\right) = \sum_{k=0}^{2} y_k\,t_k(x) = t_1 = \frac{\sin{\left(\frac{x}{2}\right)}}{\sin{\left(\frac{\pi}{4}\right)}}\frac{\sin{\left(\frac{x}{2}-\frac{\pi}{2}\right)}}{\sin{\left(-\frac{\pi}{4}\right)}} = 2 \sin{\left(\frac{x}{2}\right)}\cos{\left(\frac{x}{2}\right)} = \sin{\left(x\right)}$$

EDIT

In fact, the interpolating formula for $p\left(x\right)$, using 4 points, can be written readily as

$$p\left(x\right) = 2 s c \left[ \left(s + c \right) \frac{s c_1 - c s_1}{c_1-s_1} - \left(c-s\right) \frac{s c_3 - c s_3}{s_3 + c_3} \right]$$

with $c = \cos{\left(\frac{x}{2}\right)}$, $s = \sin{\left(\frac{x}{2}\right)}$, $c_i = \cos{\left(\frac{\alpha_i}{2}\right)}$ and $s_i = \sin{\left(\frac{\alpha_i}{2}\right)}$ as a short-hand notation.

With basic trigonometric identities $\sin{\left(x\right)} = 2 \sin{\left(\frac{x}{2}\right)} \cos{\left(\frac{x}{2}\right)}$, $\cos{\left(x\right)} = \cos^2{\left(\frac{x}{2}\right)} - \sin^2{\left(\frac{x}{2}\right)}$, this can be rewritten as

$$p\left(x\right) = \sin{\left(x\right)} \left[ 1 + \frac{1}{2} \cos{\left(x\right)} \left( \frac{s_1+c_1}{s_1-c_1} + \frac{s_3-c_3}{s_3+c_3} \right)\right]$$

So clearly, to avoid having a term $\propto \sin{\left(2 x \right)}$, the condition

$$\frac{s_1+c_1}{s_1-c_1} + \frac{s_3-c_3}{s_3+c_3} = 0$$

needs to hold. But, again, with $\alpha_1 = \frac{5}{2} \pi$ and $\alpha_3 = \frac{3}{2} \pi$, we get

$$s_1 = -\frac{\sqrt{2}}{2} \ , \ c_1 = -\frac{\sqrt{2}}{2} \\ s_3 = \frac{\sqrt{2}}{2} \ , \ c_3 = -\frac{\sqrt{2}}{2}$$

i.e. singularities in both denominators...

EDIT 2 Above condition can be rewritten to

$$\frac{\cos{\frac{1}{2}\left( \alpha_1 - \alpha_3 \right)}}{\left(\sin{\left(\frac{\alpha_1}{2}\right)}-\cos{\left(\frac{\alpha_1}{2}\right)}\right) \left(\sin{\left(\frac{\alpha_3}{2}\right)}+\cos{\left(\frac{\alpha_3}{2}\right)}\right)} = 0$$

which is satisfied for $\alpha_1-\alpha_3 = (1+2 r) \pi$.

For our simple test case, this is indeed correct, but at the same time they give infinities in the denominator, which becomes

$$\left(\sin{\left(\frac{\alpha_1}{2}\right)}-\cos{\left(\frac{\alpha_1}{2}\right)}\right)^2 = 0^2$$

But this is probably not always the case. Indeed, I have run interpolations where there were no problem. Would it just be due to an unfortunate choice of sampling points?

After some more thought I think that

1. One of the formulas on Wikipedia was slightly wrong. I have corrected it there, but the result is still $\alpha_k = \sum_{j=0, j\ne k}^{2K-1} x_j$.

2. This formula means nothing more than set the phase of the contribution of this point to the highest frequency to zero, i.e. it should be a cosine. Logically, another option is to set it to $\pi$ in order to get sines at the highest frequency. Or anything else. There is no clear preference

3. When forcing the phase to be some value, it can often happen that the denominators in the formulas contain zeros, leading to infinities. This is the case when $\sin{\frac{1}{2}\left(x_k - \alpha_k\right)} = 0$. It is very important, therefore, to have the $\alpha_k$'s satisfy the condition $x_k = \alpha_k \ne m \pi$ with $m$ an integer. This can impede the preference for the phase of the highest frequency.

4. Sometimes it is clear that the phase of the highest frequency does not matter as the amplitude should be zero, as in the case for my simple sine. There is, however, imho no criterion that can be used there, as it requires information about the function values $f_k$ and not only the values of the points $x_k$.

5. In particular, in the case for a simple sine, a phase of 0 is not possible, as this leads to infinities. An other phase has to be used, but it does not matter which: Due to the fact that $\alpha_1 -\alpha_3 = \pi$ and $y_1 = - y_3$, the global contribution to the highest frequency is zero, as required.

6. It is better just to take an odd number of points.

• Can you type the rest of your question? You can just use LaTeX syntax here for that purpose. – nicoguaro Sep 2 '16 at 21:05
• I already found it odd that I could not add equations in a simple way. Thank your for pointing that out. I will edit my question to the original format I had in mind before stack exchange started complaining about spam and I had to decrease the amount of formulas. – Toon Sep 3 '16 at 6:44
• Why not just use the nonuniform FFT (NUFFT)? It does require some oversampling but it will accomplish the goal of trigonometric interpolation you've set forth. – smh Mar 29 '19 at 12:45