Editing Trigonometric interpolation

{{Short description|Interpolation with trigonometric polynomials}}
In [[mathematics]], '''trigonometric interpolation''' is [[interpolation]] with [[trigonometric polynomial]]s. Interpolation is the process of finding a function which goes through some given [[data points]]. For trigonometric interpolation, this function has to be a trigonometric polynomial, that is, a sum of [[trigonometric function|sines and cosines]] of given periods. This form is especially suited for interpolation of [[periodic function]]s.

An important special case is when the given data points are equally spaced, in which case the solution is given by the [[discrete Fourier transform]].

==Formulation of the interpolation problem==
A trigonometric polynomial of degree ''K'' has the form
{{NumBlk|:|<math> p(x) = a_0 + \sum_{k=1}^K a_k \cos(kx) + \sum_{k=1}^K b_k \sin(kx). \,  </math>|{{EquationRef|1}}}}
This expression contains 2''K'' + 1 coefficients, ''a''<sub>0</sub>, ''a''<sub>1</sub>, … ''a''<sub>''K''</sub>, ''b''<sub>1</sub>, …, ''b''<sub>''K''</sub>, and we wish to compute those coefficients so that the function passes through ''N'' points:
:<math> p(x_n) = y_n, \quad n=0, \ldots, N-1. \, </math>
Since the trigonometric polynomial is periodic with period 2π, the ''N'' points can be distributed and ordered in one period as 
:<math> 0 \leq x_0 < x_1 < x_2 < \ldots  < x_{N-1} < 2  \pi. \, </math>
(Note that we do ''not'' in general require these points to be equally spaced.) The interpolation problem is now to find coefficients such that the trigonometric polynomial ''p'' satisfies the interpolation conditions.

==Formulation in the complex plane==
The problem becomes more natural if we formulate it in the [[complex plane]]. We can rewrite the formula for a trigonometric polynomial as
<math> p(x) = \sum_{k=-K}^K c_k e^{ikx}, \, </math>
where ''i'' is the [[imaginary unit]]. If we set ''z'' = ''e''<sup>''ix''</sup>, then this becomes
:<math> q(z) = \sum_{k=-K}^K c_k z^{k}, \, </math>
with
:<math> q(e^{ix}) \triangleq p(x). \, </math>
This reduces the problem of trigonometric interpolation to that of polynomial interpolation on the [[unit circle]]. Existence and uniqueness for trigonometric interpolation now follows immediately from the corresponding results for polynomial interpolation.

For more information on formulation of trigonometric interpolating polynomials in the complex plane, see p.&nbsp;156 of [http://www.physics.arizona.edu/~restrepo/475A/Notes/sourcea.pdf Interpolation using Fourier Polynomials].

==Solution of the problem==
Under the above conditions, there exists a solution to the problem for ''any'' given set of data points {''x''<sub>''k''</sub>, ''y''<sub>''k''</sub>} as long as ''N'', the number of data points, is not larger than the number of coefficients in the polynomial, i.e., ''N''&nbsp;≤&nbsp;2''K''+1  (a solution may or may not exist if ''N''>2''K''+1 depending upon the particular set of data points). Moreover, the  interpolating polynomial is unique if and only if the number of adjustable coefficients is equal to the number of data points, i.e., ''N''&nbsp;=&nbsp;2''K''&nbsp;+&nbsp;1. In the remainder of this article, we will assume this condition to hold true.

===Odd number of points===
If the number of points ''N'' is odd, say ''N=2K+1'', applying the [[Lagrange polynomial|Lagrange formula for polynomial interpolation]] to the polynomial formulation in the complex plane yields that the solution can be written in the form
{{NumBlk|:|<math> p(x) = \sum_{k=0}^{2K} y_k\,t_k(x),</math>|{{EquationRef|5}}}}
where
:<math> t_k(x) = e^{-iKx+iKx_k} \prod_{\begin{align}m&=0 \\[-4mu] m &\ne k\end{align}}^{2K} \frac{e^{ix}-e^{ix_m}}{e^{ix_k}-e^{ix_m}}.</math>
The factor  <math>e^{-iKx+iKx_k}</math> in this formula compensates for the fact that the complex plane formulation contains also negative powers of <math>e^{ix}</math> and is therefore not a polynomial expression in <math>e^{ix}</math>. The correctness of this expression can easily be verified by observing that <math>t_k(x_k)=1</math> and that <math>t_k(x)</math> is a linear combination of the right powers of <math>e^{ix}</math>.
Upon using the identity
{{NumBlk|:|<math>e^{iz_1}-e^{iz_2}=2i\sin\tfrac12(z_1-z_2)\,e^{(z_1 + z_2)i/2},</math>|{{EquationRef|2}}}}
the coefficient <math>t_k(x)</math> can be written in the form
{{NumBlk|:|<math> t_k(x) = \prod_{\begin{align}m&=0 \\[-4mu] m &\ne k\end{align}}^{2K} \frac{\sin\tfrac12(x-x_m)}{\sin\tfrac12(x_k-x_m)}.</math>|{{EquationRef|4}}}}

===Even number of points===
If the number of points ''N'' is even, say ''N=2K'', applying the [[Lagrange polynomial|Lagrange formula for polynomial interpolation]] to the polynomial formulation in the complex plane yields that the solution can be written in the form
{{NumBlk|:|<math> p(x) = \sum_{k=0}^{2K-1} y_k\,t_k(x),</math>|{{EquationRef|6}}}}
where
{{NumBlk|:|<math> t_k(x) = e^{-iKx+iKx_k} \frac{e^{ix}-e^{i\alpha_k}}{e^{ix_k}-e^{i\alpha_k}} \prod_{\begin{align}m&=0 \\[-4mu] m &\ne k\end{align}}^{2K-1} \frac{e^{ix}-e^{ix_m}}{e^{ix_k}-e^{ix_m}}.</math>|{{EquationRef|3}}}}
Here, the constants <math>\alpha_k</math> can be chosen freely. This is caused by the fact that the interpolating function ({{EquationNote|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 <math>\cos(Kx)</math>, i.e. the <math>\sin(Kx)</math> term vanishes, but in general the phase of the highest frequency can be chosen to be <math>\varphi_K</math>. To get an expression for <math>\alpha_k</math>, we obtain by using ({{EquationNote|2}}) that ({{EquationNote|3}}) can be written on the form
:<math> t_k(x) = \frac{\cos\tfrac12\Biggl(2Kx-\alpha_k+\displaystyle\sum\limits_{m=0,\,m \ne k}^{2K-1} x_m\Biggr)+\sum\limits_{m=-(K-1)}^{K-1}c_k e^{imx}}{2^N\sin\tfrac12(x_k-\alpha_k)\displaystyle\prod\limits_{m=0,\,m \ne k}^{2K-1}\sin\tfrac12(x_k-x_m)}.</math>
This yields
:<math>\alpha_k=\sum_{\begin{align}m&=0 \\[-4mu] m &\ne k\end{align}}^{2K-1} x_m - 2 \varphi_K</math>
and
:<math> t_k(x) = \frac{\sin\tfrac12(x-\alpha_k)}{\sin\tfrac12(x_k-\alpha_k)}\prod_{\begin{align}m&=0 \\[-4mu] m &\ne k\end{align}}^{2K-1} \frac{\sin\tfrac12(x-x_m)}{\sin\tfrac12(x_k-x_m)}.</math>

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

==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).

==Applications in numerical computing==
[[Chebfun]], a fully integrated software system written in MATLAB for computing with functions, uses trigonometric interpolation and Fourier expansions for computing with periodic functions. Many algorithms related to trigonometric interpolation are readily available in [[Chebfun]]; several examples are available [http://www.chebfun.org/docs/guide/guide11.html here].

==References==
* Kendall E. Atkinson, ''An Introduction to Numerical Analysis'' (2nd edition), Section 3.8. John Wiley & Sons, New York, 1988. {{ISBN|0-471-50023-2}}.
* M. T. Heideman, D. H. Johnson, and C. S. Burrus, "[https://web.archive.org/web/20170914220621/https://pdfs.semanticscholar.org/1790/fe007bc1ab161a1ea814748b42e3acbdc958.pdf Gauss and the history of the fast Fourier transform]," ''IEEE ASSP Magazine'' '''1''' (4), 14&ndash;21 (1984).
* G.B. Wright, M. Javed, H. Montanelli, and L.N. Trefethen, "[https://people.maths.ox.ac.uk/trefethen/trigpaper.pdf Extension of Chebfun to periodic functions]," ''SIAM. J. Sci. Comput.'', 37 (2015), C554-C573
* A. Zygmund, ''[[Trigonometric Series]]'', Volume II, Chapter X, Cambridge University Press, 1988.

==External links==
* [http://www.chebfun.org www.chebfun.org]

[[Category:Interpolation]]
[[Category:Trigonometry]]
[[Category:Articles with example MATLAB/Octave code]]