trigonometric interpolation
   HOME

TheInfoList



OR:

In
mathematics Mathematics is an area of knowledge that includes the topics of numbers, formulas and related structures, shapes and the spaces in which they are contained, and quantities and their changes. These topics are represented in modern 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 a n ...
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 at ...
. For trigonometric interpolation, this function has to be a trigonometric polynomial, that is, a sum of
sines and cosines Sines () is a city and a municipality in Portugal. The municipality, divided into two parishes, has around 14,214 inhabitants (2021) in an area of . Sines holds an important oil refinery and several petrochemical industries. It is also a popular ...
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 desc ...
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 samples of a function into a same-length sequence of equally-spaced samples of the discrete-time Fourier transform (DTFT), which is a complex- ...
.


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: : p(x_n) = y_n, \quad n=0, \ldots, N-1. \, Since the trigonometric polynomial is periodic with period 2π, the ''N'' points can be distributed and ordered in one period as : 0 \leq x_0 < x_1 < x_2 < \ldots < x_ < 2 \pi. \, (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 the ...
. We can rewrite the formula for a trigonometric polynomial as p(x) = \sum_^K c_k e^, \, 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 an ...
. If we set ''z'' = ''e''''ix'', then this becomes : q(z) = \sum_^K c_k z^, \, with : q(e^) \triangleq p(x). \, 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 Eucl ...
. 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 : t_k(x) = e^ \prod_^ \frac. The factor e^ in this formula compensates for the fact that the complex plane formulation contains also negative powers of e^ and is therefore not a polynomial expression in e^. The correctness of this expression can easily be verified by observing that t_k(x_k)=1 and that t_k(x) is a linear combination of the right powers of e^. Upon using the identity the coefficient t_k(x) 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 \alpha_k 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 \cos(Kx), i.e. the \sin(Kx) term vanishes, but in general the phase of the highest frequency can be chosen to be \varphi_K. To get an expression for \alpha_k, we obtain by using () that () can be written on the form : t_k(x) = \frac. This yields :\alpha_k=\sum_^ x_m - 2 \varphi_K and : t_k(x) = \frac\prod_^ \frac. 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 x_m are equidistant, i.e. :x_m=\frac, 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 nonneg ...
:D(x,N)=\frac +\frac \sum_^\cos(kx) = \frac, where N>0 is odd. It can easily be seen that D(x,N) is a linear combination of the right powers of e^ and satisfies :D(x_m,N)=\begin0\text m\neq0 \\1\text m=0\end. Since these two properties uniquely define the coefficients t_k(x) in (), it follows that :\begin t_k(x) &= D(x-x_k,N)=\begin \dfrac \text x\neq x_k\\ 0mu\lim\limits_ \dfrac=1 \text x= x_k \end\\&= \frac. \end 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 u ...
-function prevents any singularities and is defined by : \mathrm\,x=\frac.


Even number of points

For N even, we define the Dirichlet kernel as :D(x,N)=\frac +\frac\cos \tfrac12 Nx + \frac \sum_^\cos(kx) = \frac. Again, it can easily be seen that D(x,N) is a linear combination of the right powers of e^, does not contain the term \sin \tfrac12 Nx and satisfies :D(x_m,N)=\begin0\text m\neq0 \\1\text m=0\end. Using these properties, it follows that the coefficients t_k(x) in () are given by :\begin t_k(x) &= D(x-x_k,N)=\begin \dfrac\text x\neq x_k\\ 0mu\lim\limits_ \dfrac=1 \text x= x_k. \end\\&= \frac\cos\tfrac12 (x-x_k) \end Note that t_k(x) does not contain the \sin \tfrac12 Nx as well. Finally, note that the function \sin \tfrac12 Nx vanishes at all the points x_m. 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
here
and 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(x

0) = 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 : x_n = 2 \pi \frac, \qquad 0 \leq n < N. 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 samples of a function into a same-length sequence of equally-spaced samples of the discrete-time Fourier transform (DTFT), which is a complex- ...
(DFT) of order N. : Y_k = \sum_^ y_n \ e^ \, : y_n = p(x_n) = \frac \sum_^ Y_k \ e^ \, (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), it ...
.) 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 Lagrangiadiscrete sine transform. 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 refer ...
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 th ...
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 as a p ...
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 you ...
s,
asteroid An asteroid is a minor planet of the 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 or icy bodies with no atmosphere. ...
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/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 analogous commands for functions an ...
, 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/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 analogous commands for functions an ...
; 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