In
mathematics, trigonometric interpolation is
interpolation
In the mathematical field of numerical analysis, interpolation is a type of estimation, a method of constructing (finding) new data points based on the range of a discrete set of known data points.
In engineering and science, one often has ...
with
trigonometric polynomial In the mathematical subfields of numerical analysis and mathematical analysis, a trigonometric polynomial is a finite linear combination of functions sin(''nx'') and cos(''nx'') with ''n'' taking on the values of one or more natural numbers. The c ...
s. Interpolation is the process of finding a function which goes through some given
data points
In statistics, a unit of observation is the unit described by the data that one analyzes. A study may treat groups as a unit of observation with a country as the unit of analysis, drawing conclusions on group characteristics from data collected a ...
. For trigonometric interpolation, this function has to be a trigonometric polynomial, that is, a sum of
sines and cosines of given periods. This form is especially suited for interpolation of
periodic function
A periodic function is a function that repeats its values at regular intervals. For example, the trigonometric functions, which repeat at intervals of 2\pi radians, are periodic functions. Periodic functions are used throughout science to d ...
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
In mathematics, the discrete Fourier transform (DFT) converts a finite sequence of equally-spaced Sampling (signal processing), samples of a function (mathematics), function into a same-length sequence of equally-spaced samples of the discre ...
.
Formulation of the interpolation problem
A trigonometric polynomial of degree ''K'' has the form
This expression contains 2''K'' + 1 coefficients, ''a''
0, ''a''
1, … ''a''
''K'', ''b''
1, …, ''b''
''K'', 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 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
In mathematics, the complex plane is the plane formed by the complex numbers, with a Cartesian coordinate system such that the -axis, called the real axis, is formed by the real numbers, and the -axis, called the imaginary axis, is formed by th ...
. We can rewrite the formula for a trigonometric polynomial as
where ''i'' is the
imaginary unit
The imaginary unit or unit imaginary number () is a solution to the quadratic equation x^2+1=0. Although there is no real number with this property, can be used to extend the real numbers to what are called complex numbers, using addition a ...
. If we set ''z'' = ''e''
''ix'', then this becomes
:
with
:
This reduces the problem of trigonometric interpolation to that of polynomial interpolation on the
unit circle
In mathematics, a unit circle is a circle of unit radius—that is, a radius of 1. Frequently, especially in trigonometry, the unit circle is the circle of radius 1 centered at the origin (0, 0) in the Cartesian coordinate system in the Eucli ...
. 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 o
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 as long as ''N'', the number of data points, is not larger than the number of coefficients in the polynomial, i.e., ''N'' ≤ 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'' = 2''K'' + 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 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
and is therefore not a polynomial expression in
. The correctness of this expression can easily be verified by observing that
and that
is a linear combination of the right powers of
.
Upon using the identity
the coefficient
can be written in the form
Even number of points
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
can be chosen freely. This is caused by the fact that the interpolating function () contains an odd number of unknown constants. A common choice is to require that the highest frequency is of the form a constant times
, i.e. the
term vanishes, but in general the phase of the highest frequency can be chosen to be
. To get an expression for
, we obtain by using () that () 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.
Equidistant nodes
Further simplification of the problem is possible if nodes
are equidistant, i.e.
:
see Zygmund for more details.
Odd number of points
Further simplification by using () would be an obvious approach, but is obviously involved. A much simpler approach is to consider the
Dirichlet kernel In mathematical analysis, the Dirichlet kernel, named after the German mathematician Peter Gustav Lejeune Dirichlet, is the collection of periodic functions defined as
D_n(x)= \sum_^n e^ = \left(1+2\sum_^n\cos(kx)\right)=\frac,
where is any nonne ...
:
where
is odd. It can easily be seen that
is a linear combination of the right powers of
and satisfies
:
Since these two properties uniquely define the coefficients
in (), it follows that
:
Here, the
sinc
In mathematics, physics and engineering, the sinc function, denoted by , has two forms, normalized and unnormalized..
In mathematics, the historical unnormalized sinc function is defined for by
\operatornamex = \frac.
Alternatively, the ...
-function prevents any singularities and is defined by
:
Even number of points
For
even, we define the Dirichlet kernel as
:
Again, it can easily be seen that
is a linear combination of the right powers of
, does not contain the term
and satisfies
:
Using these properties, it follows that the coefficients
in () are given by
:
Note that
does not contain the
as well. Finally, note that the function
vanishes at all the points
. Multiples of this term can, therefore, always be added, but it is commonly left out.
Implementation
A MATLAB implementation of the above can be foun
hereand 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(x0) = 1; % fix value at x=0
Relation with the discrete Fourier transform
The special case in which the points ''x''
''n'' are equally spaced is especially important. In this case, we have
:
The transformation that maps the data points ''y''
''n'' to the coefficients ''a''
''k'', ''b''
''k'' is obtained from the
discrete Fourier transform
In mathematics, the discrete Fourier transform (DFT) converts a finite sequence of equally-spaced Sampling (signal processing), samples of a function (mathematics), function into a same-length sequence of equally-spaced samples of the discre ...
(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
In signal processing, the Nyquist frequency (or folding frequency), named after Harry Nyquist, is a characteristic of a sampler, which converts a continuous function or signal into a discrete sequence. In units of cycles per second ( Hz), i ...
.)
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
Alexis Claude Clairaut (; 13 May 1713 – 17 May 1765) was a French mathematician, astronomer, and geophysicist. He was a prominent Newtonian whose work helped to establish the validity of the principles and results that Sir Isaac Newton had out ...
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
Joseph-Louis Lagrange (born Giuseppe Luigi Lagrangia[discrete sine transform In mathematics, the discrete sine transform (DST) is a Fourier-related transform similar to the discrete Fourier transform (DFT), but using a purely real matrix. It is equivalent to the imaginary parts of a DFT of roughly twice the length, operati ...](_blank)
. The full cosine and sine interpolating polynomial, which gives rise to the DFT, was solved by
Carl Friedrich Gauss
Johann Carl Friedrich Gauss (; german: Gauß ; la, Carolus Fridericus Gauss; 30 April 177723 February 1855) was a German mathematician and physicist who made significant contributions to many fields in mathematics and science. Sometimes refe ...
in unpublished work around 1805, at which point he also derived a
fast Fourier transform
A fast Fourier transform (FFT) is an algorithm that computes the discrete Fourier transform (DFT) of a sequence, or its inverse (IDFT). Fourier analysis converts a signal from its original domain (often time or space) to a representation in t ...
algorithm to evaluate it rapidly. Clairaut, Lagrange, and Gauss were all concerned with studying the problem of inferring the
orbit
In celestial mechanics, an orbit is the curved trajectory of an object such as the trajectory of a planet around a star, or of a natural satellite around a planet, or of an artificial satellite around an object or position in space such a ...
of
planet
A planet is a large, rounded astronomical body that is neither a star nor its remnant. The best available theory of planet formation is the nebular hypothesis, which posits that an interstellar cloud collapses out of a nebula to create a ...
s,
asteroid
An asteroid is a minor planet of the Solar System#Inner solar system, inner Solar System. Sizes and shapes of asteroids vary significantly, ranging from 1-meter rocks to a dwarf planet almost 1000 km in diameter; they are rocky, metallic o ...
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
Chebfun is a Free and open source software, free/open-source software system written in MATLAB for numerical computation with functions of a real variable. It is based on the idea of overloading MATLAB's commands for vectors and matrices to analo ...
, 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
Chebfun is a Free and open source software, free/open-source software system written in MATLAB for numerical computation with functions of a real variable. It is based on the idea of overloading MATLAB's commands for vectors and matrices to analo ...
; several examples are availabl
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,
Gauss and the history of the fast Fourier transform" ''IEEE ASSP Magazine'' 1 (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
In mathematics, a trigonometric series is a infinite series of the form
: \frac+\displaystyle\sum_^(A_ \cos + B_ \sin),
an infinite version of a trigonometric polynomial.
It is called the Fourier series of the integrable function f if the term ...
'', Volume II, Chapter X, Cambridge University Press, 1988.
External links
www.chebfun.org
Interpolation
Trigonometry
Articles with example MATLAB/Octave code