Implementation A MATLAB implementation of the above can be found 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 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. 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