Skip to main content

Section 25.4 Trigonometric Interpolation

Discrete Fourier Transform allows us to interpolate any periodic function by a trigonometric polynomial by first discretizing it as \(z_k=f(k/n)\text{,}\) \(k=0, \dots, n-1\text{,}\) applying FFT to compute Fourier coefficients \(c_k\text{,}\) and then constructing a trigonometric polynomial with these coefficients. By literally following the formula (25.3.3) we arrive at the trigonometric polynomial

\begin{equation} p(x) = \frac{1}{n}\sum_{k=0}^{n-1} c_k \exp(2\pi i k x)\label{eq-bad-trig-interpolation}\tag{25.4.1} \end{equation}

This indeed interpolates the function but not in a good way, as the following code shows. It interpolates the function in Example 25.2.1 using n=9 points, which are marked as black asterisks.

f = @(x) sqrt(1 + cos(2*pi*x)) + 1i*log(2 + sin(2*pi*x));
n = 9;
z = f((0:n-1)/n);
c = fft(z);
k = 0:n-1;
p = @(t) c/n*exp(2*pi*1i*k'*t);

t = linspace(0, 1, 1000);
hold on
plot(f(t))
plot(p(t), 'r')
plot(z, 'k*')
hold off

The problem is that our formula for p involves complex exponential functions with “frequencies” \(k=0, 1, \dots, 9\) which is not how trigonometric polynomials (25.1.3) normally look. The high frequencies create the oscillatory behavior which, while matching the data, does not create a reasonable interpolant.

A key term here is aliasing, which is the fact that different trigonometric (or complex exponential) functions may coincide when restricted to a discrete grid of points (illustration). Specifically,

\begin{equation*} \exp(2\pi i k t) = \exp(2\pi i (k-n) t), \quad t = 0, 1/n, 2/n, \dots, 1 \end{equation*}

In order to avoid extraneous oscillations we should subtract \(n\) from large elements of the vector k = 0:n-1 above, so that instead of

k = [0 1 2 3 4 5 6 7 8];

we have

k = [0 1 2 3 4 -4 -3 -2 -1];

Replacing k = 0:n-1 in the above code with k = [0:(n-1)/2 (1-n)/2:-1]; corrects the issue.

The previous paragraph operates with an odd number \(n\) which is more natural in this context because trigonometric polynomials (25.1.3) normally contain an odd number of terms. If we have an even number of points to interpolate through, the aliasing gets more involved: for example, with \(n=6\) points we have frequencies 0, 1, 2, 3, 4, 5 and it is not clear whether to keep \(k=3\) as is or to replace it with \(k-n=-3\text{.}\) The best approach is to split this term in two, so that the list of coefficients \(k\) becomes 0, 1, 2, 3, -3, -2, -1. The corresponding Fourier coefficient is split in two halves accordingly:

c = [c(1:n/2) c(n/2+1)/2 c(n/2+1)/2 c(n/2+2:end)];    
k = [0:n/2 -n/2:-1];

If it is not known whether \(n\) is even or odd, one can use a conditional statement to handle both cases:

if mod(n, 2) == 1
    k = [0:(n-1)/2 (1-n)/2:-1];
else
    c = [c(1:n/2) c(n/2+1)/2 c(n/2+1)/2 c(n/2+2:end)];
    k = [0:n/2 -n/2:-1];
end

Choose a few points in the complex plane so that trigonometric interpolation through these points creates a curve similar to Meta logo.

Answer

It is natural to place the self-intersection point at \(0\text{,}\) then the points of interpolation will be symmetric about the vertical (imaginary) axis. Trial and error leads to the choice

z = [6, 6-6i, 2-3i, -2+3i, -6, -6-6i, -2-3i, 2+3i];

which describes the curve starting from the upper right part, going down, then diagonally to top left, then down, and diagonally to top right again. We have \(n=8\) points, so the above discussion about odd/even distinction is relevant. For the sake of reusability, the code uses a conditional statement that can handle both odd and even cases.

n = numel(z);
c = fft(z);
if mod(n, 2) == 1
    k = [0:(n-1)/2 (1-n)/2:-1]; 
else
    c = [c(1:n/2) c(n/2+1)/2 c(n/2+1)/2 c(n/2+2:end)];
    k = [0:n/2 -n/2:-1];
end
p = @(t) c/n*exp(2*pi*1i*k'*t);

t = linspace(0, 1, 1000);
plot(p(t), 'linewidth', 30)
hold on
plot(z, 'k.', 'markersize', 30)
hold off
axis equal

The curve is made extra thick because this is how Meta logo looks. The code also marks the eight points of interpolation.

The preceding examples had complex-valued data, but the same approach allows us to interpolate real-valued data with a real trigonometric polynomial. The main difference is in plotting the real-valued periodic function: we want the horizontal axis to represent \(t\) while the vertical axis shows the values of the function.

Plot a trigonometric polynomial with period \(1\) which attains prescribed (random) real values at the points \(0, 0.2, 0.4, 0.6, 0.8\text{.}\)

Answer

Here \(n=5\) (the number of interpolation points), and z = rand(1, n) (real numbers). Theoretically, the interpolating polynomial will be real-valued but due to limited precision, some small imaginary component may appear. So we should take the real part real(p(t)).

n = 5;
z = rand(1, n);
c = fft(z);
k = [0:(n-1)/2 (1-n)/2:-1];
p = @(t) c/n*exp(2*pi*1i*k'*t);

t = linspace(0, 1, 1000);
plot(t, real(p(t)), (0:n-1)/n, z, 'r*')

Note that \(p(1)=p(0)\) by construction.