Editing Trigonometric interpolation (section)

==Equidistant nodes==
Further simplification of the problem is possible if nodes <math>x_m</math> are equidistant, i.e.
:<math>x_m=\frac{2\pi m}{N},</math>
see Zygmund for more details.

===Odd number of points===
Further simplification by using ({{EquationNote|4}}) would be an obvious approach, but is obviously involved. A much simpler approach is to consider the [[Dirichlet kernel]]
:<math>D(x,N)=\frac{1}{N} +\frac{2}{N} \sum_{k=1}^{(N-1)/2}\cos(kx) = \frac{\sin\tfrac12 Nx}{N\sin\tfrac12 x},</math>
where <math>N>0</math> is odd. It can easily be seen that <math>D(x,N)</math> is a linear combination of the right powers of <math>e^{ix}</math> and satisfies
:<math>D(x_m,N)=\begin{cases}0\text{ for } m\neq0 \\1\text{ for } m=0\end{cases}.</math>
Since these two properties uniquely define the coefficients <math>t_k(x)</math> in ({{EquationNote|5}}), it follows that
:<math>\begin{align}
t_k(x) &= D(x-x_k,N)=\begin{cases}
\dfrac{\sin\tfrac12 N(x-x_k)}{N\sin\tfrac12 (x-x_k)} \text{ for } x\neq x_k\\[10mu]
\lim\limits_{x\to 0} \dfrac{\sin\tfrac12 Nx}{N\sin\tfrac12 x}=1 \text{ for } x= x_k
\end{cases}\\&= \frac{\mathrm{sinc}\,\tfrac12 N(x-x_k)}{\mathrm{sinc}\,\tfrac12 (x-x_k)}.
\end{align}</math>
Here, the [[sinc]]-function prevents any singularities and is defined by
:<math> \mathrm{sinc}\,x=\frac{\sin x}{x}.</math>

===Even number of points===
For <math>N</math> even, we define the Dirichlet kernel as
:<math>D(x,N)=\frac{1}{N} +\frac{1}{N}\cos \tfrac12 Nx + \frac{2}{N} \sum_{k=1}^{(N-1)/2}\cos(kx) = \frac{\sin\tfrac12 Nx}{N\tan\tfrac12 x}.</math>
Again, it can easily be seen that <math>D(x,N)</math> is a linear combination of the right powers of <math>e^{ix}</math>, does not contain the term <math> \sin \tfrac12 Nx </math> and satisfies
:<math>D(x_m,N)=\begin{cases}0\text{ for } m\neq0 \\1\text{ for } m=0\end{cases}.</math>
Using these properties, it follows that the coefficients <math>t_k(x)</math> in ({{EquationNote|6}}) are given by
:<math>\begin{align}
t_k(x) &= D(x-x_k,N)=\begin{cases}
\dfrac{\sin\tfrac12 N(x-x_k)}{N\tan\tfrac12 (x-x_k)}\text{ for } x\neq x_k\\[10mu]
\lim\limits_{x\to 0} \dfrac{\sin\tfrac12 Nx}{N\tan\tfrac12 x}=1 \text{ for } x= x_k.
\end{cases}\\&= \frac{\mathrm{sinc}\,\tfrac12 N(x-x_k)}{ \mathrm{sinc}\,\tfrac12 (x-x_k)}\cos\tfrac12 (x-x_k)
\end{align}</math>
Note that <math>t_k(x)</math> does not contain the <math> \sin \tfrac12 Nx </math> as well. Finally, note that the function <math> \sin \tfrac12 Nx </math> vanishes at all the points <math>x_m</math>. Multiples of this term can, therefore, always be added, but it is commonly left out.

===Implementation===
A MATLAB implementation of the above can be found [http://www.math.udel.edu/~braun/M428/Matlab/interpolation/triginterp.m here] and is given by:
<syntaxhighlight lang="matlab">
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
</syntaxhighlight>

===Relation with the discrete Fourier transform===
The special case in which the points ''x''<sub>''n''</sub> are equally spaced is especially important. In this case, we have
:<math> x_n = 2 \pi \frac{n}{N}, \qquad 0 \leq n < N.</math>

The transformation that maps the data points ''y''<sub>''n''</sub> to the coefficients ''a''<sub>''k''</sub>, ''b''<sub>''k''</sub> is obtained from the [[discrete Fourier transform]] (DFT) of order&nbsp;N.

:<math> Y_k = \sum_{n=0}^{N-1} y_n \ e^{-i 2 \pi nk/N} \, </math>

:<math> y_n = p(x_n) = \frac{1}{N} \sum_{k=0}^{N-1} Y_k \ e^{i 2 \pi nk/N} \, </math>

(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 and odd functions|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 [[Cooley–Tukey FFT algorithm|fast Fourier transform]] algorithm to evaluate it rapidly. Clairaut, Lagrange, and Gauss were all concerned with studying the problem of inferring the [[orbit]] of [[planet]]s, [[asteroid]]s, 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).