Interpolation is the process of finding a function which goes through some given data points.
This form is especially suited for interpolation of periodic functions.
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.
A trigonometric polynomial of degree K has the form This expression contains 2K + 1 coefficients, a0, a1, … aK, b1, …, bK, and we wish to compute those coefficients so that the function passes through N points: Since the trigonometric polynomial is periodic with period 2π, the N points can be distributed and ordered in one period as (Note that we do not in general require these points to be equally spaced.)
The problem becomes more natural if we formulate it in the complex plane.
If we set z = eix, then this becomes with This reduces the problem of trigonometric interpolation to that of polynomial interpolation on the unit circle.
For more information on formulation of trigonometric interpolating polynomials in the complex plane, see p. 156 of Interpolation using Fourier Polynomials.
Under the above conditions, there exists a solution to the problem for any given set of data points {xk, yk} as long as N, 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 if N>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 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 where The factor
in this formula compensates for the fact that the complex plane formulation contains also negative powers of
The correctness of this expression can easily be verified by observing that
can be written in the form 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 where Here, the constants
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
term vanishes, but in general the phase of the highest frequency can be chosen to be
, we obtain by using (2) that (3) can be written on the form This yields and Note that care must be taken in order to avoid infinities caused by zeros in the denominators.
and satisfies Since these two properties uniquely define the coefficients
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: The special case in which the points xn are equally spaced is especially important.
In this case, we have The transformation that maps the data points yn to the coefficients ak, bk is obtained from the discrete Fourier transform (DFT) of order 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.
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 here.