Trigonometric interpolation

If the number of points N is odd, say N=2K+1, applying the Lagrange formula for polynomial interpolation to the polynomial formulation in the complex plane yields that the solution can be written in the form

${\displaystyle p(x)=\sum _{k=0}^{2K}y_{k}\,t_{k}(x),}$

(5)

where

${\displaystyle t_{k}(x)=e^{-iKx+iKx_{k}}\prod _{\begin{aligned}m&=0\\[-4mu]m&\neq k\end{aligned}}^{2K}{\frac {e^{ix}-e^{ix_{m}}}{e^{ix_{k}}-e^{ix_{m}}}}.}$

The factor ${\displaystyle e^{-iKx+iKx_{k}}}$ in this formula compensates for the fact that the complex plane formulation contains also negative powers of ${\displaystyle e^{ix}}$ and is therefore not a polynomial expression in ${\displaystyle e^{ix}}$. The correctness of this expression can easily be verified by observing that ${\displaystyle t_{k}(x_{k})=1}$ and that ${\displaystyle t_{k}(x)}$ is a linear combination of the right powers of ${\displaystyle e^{ix}}$. Upon using the identity

${\displaystyle e^{iz_{1}}-e^{iz_{2}}=2i\sin {\tfrac {1}{2}}(z_{1}-z_{2})\,e^{(z_{1}+z_{2})i/2},}$

(2)

the coefficient ${\displaystyle t_{k}(x)}$ can be written in the form

${\displaystyle t_{k}(x)=\prod _{\begin{aligned}m&=0\\[-4mu]m&\neq k\end{aligned}}^{2K}{\frac {\sin {\tfrac {1}{2}}(x-x_{m})}{\sin {\tfrac {1}{2}}(x_{k}-x_{m})}}.}$

(4)

If the number of points N is even, say N=2K, applying the Lagrange formula for polynomial interpolation to the polynomial formulation in the complex plane yields that the solution can be written in the form

${\displaystyle p(x)=\sum _{k=0}^{2K-1}y_{k}\,t_{k}(x),}$

(6)

where

${\displaystyle t_{k}(x)=e^{-iKx+iKx_{k}}{\frac {e^{ix}-e^{i\alpha _{k}}}{e^{ix_{k}}-e^{i\alpha _{k}}}}\prod _{\begin{aligned}m&=0\\[-4mu]m&\neq k\end{aligned}}^{2K-1}{\frac {e^{ix}-e^{ix_{m}}}{e^{ix_{k}}-e^{ix_{m}}}}.}$

(3)

Here, the constants ${\displaystyle \alpha _{k}}$ can be chosen freely. This is caused by the fact that the interpolating function (1) contains an odd number of unknown constants. A common choice is to require that the highest frequency is of the form a constant times ${\displaystyle \cos(Kx)}$, i.e. the ${\displaystyle \sin(Kx)}$ term vanishes, but in general the phase of the highest frequency can be chosen to be ${\displaystyle \varphi _{K}}$. To get an expression for ${\displaystyle \alpha _{k}}$, we obtain by using (2) that (3) can be written on the form

${\displaystyle t_{k}(x)={\frac {\cos {\tfrac {1}{2}}{\Biggl (}2Kx-\alpha _{k}+\displaystyle \sum \limits _{m=0,\,m\neq k}^{2K-1}x_{m}{\Biggr )}+\sum \limits _{m=-(K-1)}^{K-1}c_{k}e^{imx}}{2^{N}\sin {\tfrac {1}{2}}(x_{k}-\alpha _{k})\displaystyle \prod \limits _{m=0,\,m\neq k}^{2K-1}\sin {\tfrac {1}{2}}(x_{k}-x_{m})}}.}$

This yields

${\displaystyle \alpha _{k}=\sum _{\begin{aligned}m&=0\\[-4mu]m&\neq k\end{aligned}}^{2K-1}x_{m}-2\varphi _{K}}$

and

${\displaystyle t_{k}(x)={\frac {\sin {\tfrac {1}{2}}(x-\alpha _{k})}{\sin {\tfrac {1}{2}}(x_{k}-\alpha _{k})}}\prod _{\begin{aligned}m&=0\\[-4mu]m&\neq k\end{aligned}}^{2K-1}{\frac {\sin {\tfrac {1}{2}}(x-x_{m})}{\sin {\tfrac {1}{2}}(x_{k}-x_{m})}}.}$

Note that care must be taken in order to avoid infinities caused by zeros in the denominators.

Further simplification by using (4) would be an obvious approach, but is obviously involved. A much simpler approach is to consider the Dirichlet kernel

${\displaystyle D(x,N)={\frac {1}{N}}+{\frac {2}{N}}\sum _{k=1}^{(N-1)/2}\cos(kx)={\frac {\sin {\tfrac {1}{2}}Nx}{N\sin {\tfrac {1}{2}}x}},}$

where ${\displaystyle N>0}$ is odd. It can easily be seen that ${\displaystyle D(x,N)}$ is a linear combination of the right powers of ${\displaystyle e^{ix}}$ and satisfies

${\displaystyle D(x_{m},N)={\begin{cases}0{\text{ for }}m\neq 0\\1{\text{ for }}m=0\end{cases}}.}$

Since these two properties uniquely define the coefficients ${\displaystyle t_{k}(x)}$ in (5), it follows that

${\displaystyle {\begin{aligned}t_{k}(x)&=D(x-x_{k},N)={\begin{cases}{\dfrac {\sin {\tfrac {1}{2}}N(x-x_{k})}{N\sin {\tfrac {1}{2}}(x-x_{k})}}{\text{ for }}x\neq x_{k}\\[10mu]\lim \limits _{x\to 0}{\dfrac {\sin {\tfrac {1}{2}}Nx}{N\sin {\tfrac {1}{2}}x}}=1{\text{ for }}x=x_{k}\end{cases}}\\&={\frac {\mathrm {sinc} \,{\tfrac {1}{2}}N(x-x_{k})}{\mathrm {sinc} \,{\tfrac {1}{2}}(x-x_{k})}}.\end{aligned}}}$

Here, the sinc-function prevents any singularities and is defined by

${\displaystyle \mathrm {sinc} \,x={\frac {\sin x}{x}}.}$

For ${\displaystyle N}$ even, we define the Dirichlet kernel as

${\displaystyle D(x,N)={\frac {1}{N}}+{\frac {1}{N}}\cos {\tfrac {1}{2}}Nx+{\frac {2}{N}}\sum _{k=1}^{(N-1)/2}\cos(kx)={\frac {\sin {\tfrac {1}{2}}Nx}{N\tan {\tfrac {1}{2}}x}}.}$

Again, it can easily be seen that ${\displaystyle D(x,N)}$ is a linear combination of the right powers of ${\displaystyle e^{ix}}$, does not contain the term ${\displaystyle \sin {\tfrac {1}{2}}Nx}$ and satisfies

${\displaystyle D(x_{m},N)={\begin{cases}0{\text{ for }}m\neq 0\\1{\text{ for }}m=0\end{cases}}.}$

Using these properties, it follows that the coefficients ${\displaystyle t_{k}(x)}$ in (6) are given by

${\displaystyle {\begin{aligned}t_{k}(x)&=D(x-x_{k},N)={\begin{cases}{\dfrac {\sin {\tfrac {1}{2}}N(x-x_{k})}{N\tan {\tfrac {1}{2}}(x-x_{k})}}{\text{ for }}x\neq x_{k}\\[10mu]\lim \limits _{x\to 0}{\dfrac {\sin {\tfrac {1}{2}}Nx}{N\tan {\tfrac {1}{2}}x}}=1{\text{ for }}x=x_{k}.\end{cases}}\\&={\frac {\mathrm {sinc} \,{\tfrac {1}{2}}N(x-x_{k})}{\mathrm {sinc} \,{\tfrac {1}{2}}(x-x_{k})}}\cos {\tfrac {1}{2}}(x-x_{k})\end{aligned}}}$

Note that ${\displaystyle t_{k}(x)}$ does not contain the ${\displaystyle \sin {\tfrac {1}{2}}Nx}$ as well. Finally, note that the function ${\displaystyle \sin {\tfrac {1}{2}}Nx}$ vanishes at all the points ${\displaystyle x_{m}}$. Multiples of this term can, therefore, always be added, but it is commonly left out.

A MATLAB implementation of the above can be found here and is given by:

function P = triginterp(xi,x,y)
% TRIGINTERP Trigonometric interpolation.
% Input:
%   xi  evaluation points for the interpolant (vector)
%   x   equispaced interpolation nodes (vector, length N)
%   y   interpolation values (vector, length N)
% Output:
%   P   values of the trigonometric interpolant (vector)
N = length(x);
% Adjust the spacing of the given independent variable.
h = 2/N;
scale = (x(2)-x(1)) / h;
x = x/scale;  xi = xi/scale;
% Evaluate interpolant.
P = zeros(size(xi));
for k = 1:N
  P = P + y(k)*trigcardinal(xi-x(k),N);
end

function tau = trigcardinal(x,N)
ws = warning('off','MATLAB:divideByZero');
% Form is different for even and odd N.
if rem(N,2)==1   % odd
  tau = sin(N*pi*x/2) ./ (N*sin(pi*x/2));
else             % even
  tau = sin(N*pi*x/2) ./ (N*tan(pi*x/2));
end
warning(ws)
tau(x==0) = 1;     % fix value at x=0

The special case in which the points x_n are equally spaced is especially important. In this case, we have

${\displaystyle x_{n}=2\pi {\frac {n}{N}},\qquad 0\leq n<N.}$

The transformation that maps the data points y_n to the coefficients a_k, b_k is obtained from the discrete Fourier transform (DFT) of order N.

${\displaystyle Y_{k}=\sum _{n=0}^{N-1}y_{n}\ e^{-i2\pi nk/N}\,}$

${\displaystyle y_{n}=p(x_{n})={\frac {1}{N}}\sum _{k=0}^{N-1}Y_{k}\ e^{i2\pi nk/N}\,}$

(Because of the way the problem was formulated above, we have restricted ourselves to odd numbers of points. This is not strictly necessary; for even numbers of points, one includes another cosine term corresponding to the Nyquist frequency.)

The case of the cosine-only interpolation for equally spaced points, corresponding to a trigonometric interpolation when the points have even symmetry, was treated by Alexis Clairaut in 1754. In this case the solution is equivalent to a discrete cosine transform. The sine-only expansion for equally spaced points, corresponding to odd symmetry, was solved by Joseph Louis Lagrange in 1762, for which the solution is a discrete sine transform. The full cosine and sine interpolating polynomial, which gives rise to the DFT, was solved by Carl Friedrich Gauss in unpublished work around 1805, at which point he also derived a fast Fourier transform algorithm to evaluate it rapidly. Clairaut, Lagrange, and Gauss were all concerned with studying the problem of inferring the orbit of planets, asteroids, etc., from a finite set of observation points; since the orbits are periodic, a trigonometric interpolation was a natural choice. See also Heideman et al. (1984).