Inmathematics,trigonometric interpolation isinterpolation withtrigonometric polynomials. Interpolation is the process of finding a function which goes through some givendata points. For trigonometric interpolation, this function has to be a trigonometric polynomial, that is, a sum ofsines and cosines of given periods. This form is especially suited for interpolation ofperiodic functions.

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

A trigonometric polynomial of degreeK has the form

${\displaystyle p(x)=a_0+\sum _{k=1}^K{a}_k\cos (kx)+\sum _{k=1}^K{b}_k\sin (kx).}$

1

This expression contains 2K + 1 coefficients,a₀,a₁, …a_K,b₁, …,b_K, and we wish to compute those coefficients so that the function passes throughN points:

${\displaystyle p(x_n)=y_n,n=0,\dots ,N-1.}$

Since the trigonometric polynomial is periodic with period 2π, theN points can be distributed and ordered in one period as

(Note that we donot in general require these points to be equally spaced.) The interpolation problem is now to find coefficients such that the trigonometric polynomialp satisfies the interpolation conditions.

The problem becomes more natural if we formulate it in thecomplex plane. We can rewrite the formula for a trigonometric polynomial as${\displaystyle p(x)=\sum _{k=-K}^K{c}_k{e}^{ikx},}$wherei is theimaginary unit. If we setz =e^ix, then this becomes

${\displaystyle q(z)=\sum _{k=-K}^K{c}_k{z}^k,}$

with

${\displaystyle q(e^{ix})\triangleq p(x).}$

This reduces the problem of trigonometric interpolation to that of polynomial interpolation on theunit 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. 156 ofInterpolation using Fourier Polynomials.

Under the above conditions, there exists a solution to the problem forany given set of data points {x_k,y_k} as long asN, the number of data points, is not larger than the number of coefficients in the polynomial, i.e.,N ≤ 2K+1 (a solution may or may not exist ifN>2K+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 = 2K + 1. In the remainder of this article, we will assume this condition to hold true.

If the number of pointsN is odd, sayN=2K+1, applying theLagrange 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

The factor${\displaystyle e^{-iKx+iK{x}_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^{i{z}_1}-e^{i{z}_2}=2i\sin \frac{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

4

If the number of pointsN is even, sayN=2K, applying theLagrange 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

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 {\phi }_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 \frac{1}{2}(2Kx-{\alpha }_k+{\displaystyle \sum _{m=0,m\ne k}^{2K-1}{x}_m)+\sum _{m=-(K-1)}^{K-1}{c}_k{e}^{imx}}}{2^N\sin \frac{1}{2}(x_k-{\alpha }_k){\displaystyle \prod _{m=0,m\ne k}^{2K-1}\sin \frac{1}{2}(x_k-x_m)}}.}$

This yields

and

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

Further simplification of the problem is possible if nodes${\displaystyle x_m}$ are equidistant, i.e.

${\displaystyle x_m=\frac{2\pi m}{N},}$

see Zygmund for more details.

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

${\displaystyle D(x,N)=\frac{1}{N}+\frac{2}{N}\sum _{k=1}^{(N-1)/2}\cos (kx)=\frac{\sin \frac{1}{2}Nx}{N\sin \frac{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{array}{l}0\text{ for }m\ne 0\\ 1\text{ for }m=0\end{array}.}$

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

Here, thesinc-function prevents any singularities and is defined by

${\displaystyle \mathrm{s}\mathrm{i}\mathrm{n}\mathrm{c}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 \frac{1}{2}Nx+\frac{2}{N}\sum _{k=1}^{(N-1)/2}\cos (kx)=\frac{\sin \frac{1}{2}Nx}{N\tan \frac{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 \frac{1}{2}Nx}$ and satisfies

${\displaystyle D(x_m,N)=\{\begin{array}{l}0\text{ for }m\ne 0\\ 1\text{ for }m=0\end{array}.}$

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

Note that${\displaystyle t_k(x)}$ does not contain the${\displaystyle \sin \frac{1}{2}Nx}$ as well. Finally, note that the function${\displaystyle \sin \frac{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.

AMATLAB implementation of the above can be foundhere and is given by:

functionP=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));fork=1:NP=P+y(k)*trigcardinal(xi-x(k),N);endfunctiontau=trigcardinal(x,N)ws=warning('off','MATLAB:divideByZero');% Form is different for even and odd N.ifrem(N,2)==1% oddtau=sin(N*pi*x/2)./(N*sin(pi*x/2));else% eventau=sin(N*pi*x/2)./(N*tan(pi*x/2));endwarning(ws)tau(x==0)=1;% fix value at x=0

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

The transformation that maps the data pointsy_n to the coefficientsa_k,b_k is obtained from thediscrete 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 theNyquist frequency.)

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

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 inChebfun; several examples are availablehere.

• 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, "Gauss and the history of the fast Fourier transform,"IEEE ASSP Magazine1 (4), 14–21 (1984).
• G.B. Wright, M. Javed, H. Montanelli, and L.N. Trefethen, "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.

• www.chebfun.org

• Interpolation
• Trigonometry

• Articles with short description
• Short description is different from Wikidata
• Articles with example MATLAB/Octave code