Chapter 3 - Trigonometric Interpolation

In the last chapter, we found a powerful tool for interpolation: the polynomial. With a careful choice of points, the Chebyshev points, we could create excellent approximations. But polynomials have a fundamental character: they want to go to infinity. They are not periodic. If you try to model a phenomenon that repeats itself, a sound wave, a heartbeat, the seasons, forcing a polynomial to do the job is like trying to draw a circle with a set of straight rulers. It can be done, but it’s awkward and unnatural.

We need a new set of building blocks, a new language designed from the ground up for periodicity. This is the world of sines and cosines. This is the world of Jean-Baptiste Joseph Fourier. His revolutionary idea, born from studying the flow of heat, was that any periodic function could be represented as a sum of simple, harmonically related sinusoids. This was more than a new technique; it was a new way of seeing. It gave us a pair of glasses to switch between viewing a function as a shape in time and as a spectrum of frequencies.

The Continuous World: Fourier Series

Instead of building our functions from powers of $x$, like $\{1,x,x^2,\dots \}$, we will build them from a basis of pure oscillations. It’s most elegant to use complex exponentials, which wrap up sine and cosine into one beautiful package thanks to Euler’s formula, $e^{i\theta }=\cos \theta +i\sin \theta$. Using $e^{i\theta }$ is mathematically convenient because it combines the properties of both sine and cosine, simplifying the algebra greatly.

A trigonometric polynomial of degree $m$ is a function of the form:

$$p_m(t)=\sum _{j=-m}^m{c}_j{e}^{2\pi ijt}$$

where the $c_j$ are complex coefficients. This function is naturally periodic with a period of 1. If the coefficients satisfy the symmetry $\overline{c_j}=c_{-j}$ (which means the coefficient for frequency $j$ is the complex conjugate of the coefficient for frequency $-j$), the function will be purely real-valued.

Functions as Vectors in an Infinite-Dimensional Space

To truly grasp Fourier’s idea, we must think like a linear algebraist, as Gilbert Strang would urge us to. Let’s consider the space of all complex-valued functions on the interval $[0,1]$ that have a finite “energy,” a space called $L^2(0,1)$. (Think of $L^2$ as the space of all functions that have a finite “length” or magnitude when measured continuously.) In this space, functions are not just curves; they are vectors. We can define an inner product (the continuous version of a dot product) between two functions $f$ and $g$:

$$⟨f,g⟩={\int }_0^{1}f(t)\overline{g(t)}dt$$

This inner product measures how much the two functions align or overlap. The “length” or norm of a function is then $\Vert f{\Vert }_{L^2}=\sqrt{⟨f,f⟩}$.

The magic of the complex exponentials, ${\varphi }_k(t)=e^{2\pi ikt}$, is that they form an orthonormal basis for this infinite-dimensional vector space. “Orthonormal” means these building blocks are perfectly independent of each other (orthogonal) and all have unit length (normalized). Let’s check the orthogonality:

$$⟨{\varphi }_k,{\varphi }_j⟩={\int }_0^1{e}^{2\pi ikt}\overline{e^{2\pi ijt}}dt={\int }_0^1{e}^{2\pi i(k-j)t}dt=\{\begin{array}{ll}1& \text{if }k=j\\ 0& \text{if }k\ne j\end{array}$$

They are perfectly orthogonal, just like the standard basis vectors ${\mathbf{e}}_1,{\mathbf{e}}_2,\dots$ in ${\mathbb{R}}^n$. This orthogonality is not just a mathematical curiosity; it is the key that unlocks everything.

Finding the Fourier Coefficients

Because the basis is orthonormal, finding the coordinates of a function is easy. Just as you find the $k$-th coordinate of a vector $\mathbf{v}$ by projecting it onto the $k$-th basis vector ($\mathbf{v}\cdot {\mathbf{e}}_k$), we find the Fourier coefficients $\hat{f}(k)$ of a function $f(t)$ by projecting it onto our basis functions:

$$\hat{f}(k)=⟨f,{\varphi }_k⟩={\int }_0^{1}f(t)e^{-2\pi ikt}dt$$

The coefficient $\hat{f}(k)$ tells us exactly how much of the wave with frequency $k$ is contained in the function $f$.

The Fourier Series Theorem then states that we can reconstruct the original function perfectly by summing up these basis functions weighted by their coefficients:

$$f(t)=\sum _{k=-\infty }^{\infty }\hat{f}(k)e^{2\pi ikt}$$

This is a profound statement. It says any well-behaved periodic function is a superposition of simple waves. The set of coefficients $\{\hat{f}(k)\}$ is the function’s representation in the frequency domain.

The Dialogue Between Smoothness and Decay

Here we come to a crucial insight. There is a beautiful and deep connection between how smooth a function is in the time domain and how quickly its Fourier coefficients decay in the frequency domain. What happens when we differentiate a Fourier series?

$$f^{\prime }(t)=\sum _{k=-\infty }^{\infty }(2\pi ik)\hat{f}(k)e^{2\pi ikt}$$

The Fourier coefficients of the derivative, $\widehat{f^{\prime }}(k)$, are simply $2\pi ik\cdot \hat{f}(k)$. For the series of the derivative to converge (i.e., for its energy to be finite), the coefficients $\hat{f}(k)$ must decay faster than $1/k$. For the second derivative, they must decay faster than $1/k^2$, and so on.

Smoothness is Rapid Decay

The smoother a function is, the faster its Fourier coefficients $\hat{f}(k)$ must go to zero as the frequency $|k|\to \infty$. This means for smooth functions, we only need to keep a few low-frequency coefficients to achieve excellent accuracy.

• A function with a jump discontinuity (the least smooth): $\hat{f}(k)\sim O(1/k)$.
• A continuous function with a sharp corner: $\hat{f}(k)\sim O(1/k^2)$.
• An infinitely smooth function (like a Gaussian): $\hat{f}(k)$ decays faster than any power of $1/k$ (exponential decay).

This dialogue is the key to signal compression, denoising, and understanding why our methods work.

The Discrete World: The Discrete Fourier Transform (DFT)

The Fourier series is a tool for the continuous world of pure mathematics. On a computer, we only have a finite number of samples. Suppose we measure a function $f(t)$ at $N$ equally spaced points in the interval $[0,1)$:

$$t_j=\frac{j}{N},j=0,1,\dots ,N-1$$

Our data is a vector of values $\mathbf{y}=(y_0,y_1,\dots ,y_{N-1})$, where $y_j=f(t_j)$.

How can we approximate the Fourier coefficients from this discrete data? The most natural approach is to replace the integral in the definition of $\hat{f}(k)$ with a numerical quadrature rule. The simplest one for a periodic function is the trapezoidal rule (which, in this case, is equivalent to a simple Riemann sum):

$$\hat{f}(k)={\int }_0^{1}f(t)e^{-2\pi ikt}dt\approx \frac{1}{N}\sum _{j=0}^{N-1}f(t_j)e^{-2\pi ik{t}_j}$$

Let’s define a vector of coefficients, $\mathbf{c}$, based on this sum (we’ll ignore the $1/N$ scaling factor for now):

$$c_k=\sum _{j=0}^{N-1}{y}_j{e}^{-2\pi ikj/N}$$

This formula is the Discrete Fourier Transform (DFT). It is a linear transformation that takes a vector of $N$ samples in the “time domain” and produces a vector of $N$ coefficients in the “frequency domain.”

The Miracle of Discrete Interpolation

Now for a moment of magic. What happens if we build a trigonometric polynomial using these approximated coefficients, $c_k/N$?

$$p(t)=\frac{1}{N}\sum _{k=0}^{N-1}{c}_k{e}^{2\pi ikt}$$

Let’s evaluate this polynomial at one of our sample points, $t_l=l/N$:

$$p(t_l)=\frac{1}{N}\sum _{k=0}^{N-1}\left(\sum _{j=0}^{N-1}{y}_j{e}^{-2\pi ikj/N}\right)e^{2\pi ikl/N}=\sum _{j=0}^{N-1}{y}_j\left(\frac{1}{N}\sum _{k=0}^{N-1}{e}^{2\pi ik(l-j)/N}\right)$$

The term in the parenthesis is a geometric sum of the powers of the $N$-th root of unity, ${\omega }_N=e^{2\pi i/N}$. This sum has a remarkable property, often called the discrete orthogonality relation:

$$\frac{1}{N}\sum _{k=0}^{N-1}({\omega }_N^{l-j}{)}^k=\{\begin{array}{ll}1& \text{if }l-j\text{ is a multiple of }N\\ 0& \text{otherwise}\end{array}$$

(Don’t worry about the detailed derivation here; the key is the outcome.) Since both $l$ and $j$ are between $0$ and $N-1$, this is equivalent to saying the sum is 1 if $l=j$ and 0 otherwise. The inner sum collapses, acting like a Kronecker delta, picking out only the term where $j=l$. The result is:

$$p(t_l)=y_l$$

This is a stunning result. The trigonometric polynomial built from the DFT coefficients exactly interpolates the original data points. The numerical approximation we made for the continuous integral has, by a beautiful symmetry of mathematics (the discrete orthogonality), become an exact interpolation formula.

The DFT as a Change of Basis

The DFT is a linear transformation, which we can write as a matrix-vector product, $\mathbf{c}=F_N\mathbf{y}$. The matrix $F_N$ is the DFT matrix:

$$(F_N{)}_{kj}={\omega }_N^{-kj}=e^{-2\pi ikj/N}$$

This matrix is a change of basis matrix. It transforms a vector $\mathbf{y}$ from the standard basis (the “time domain”) into the discrete Fourier basis (the “frequency domain”). The columns of the inverse Fourier matrix, $F_N^{-1}=\frac{1}{N}{F}_N^H$, are the discrete basis vectors themselves.

The Algorithm: The Fast Fourier Transform (FFT)

The DFT is a beautiful theoretical tool, but a direct matrix-vector multiplication (calculating all $N$ coefficients by summing $N$ terms each) costs $O(N^2)$ operations. For a million data points ($N=10^6$), that’s a trillion operations, impossibly slow.

The Fast Fourier Transform (FFT) is a revolutionary algorithm that computes the exact same DFT in only $O(N\log N)$ time. For $N=10^6$, this is the difference between a calculation taking hours or days and one taking milliseconds. It was rediscovered by Cooley and Tukey in 1965 and is one of the most important algorithms ever developed.

The idea is a classic divide and conquer. A DFT of size $N$ can be broken down into smaller DFTs. If $N=2m$ (an even number), we can split the sum into even and odd indices:

$$\begin{array}{rl}{c}_k& =\sum _{j=0}^{m-1}{y}_{2j}{e}^{-2\pi ik(2j)/N}+\sum _{j=0}^{m-1}{y}_{2j+1}{e}^{-2\pi ik(2j+1)/N}\\ & =\sum _{j=0}^{m-1}{y}_{2j}(e^{-2\pi ikj/m})+e^{-2\pi ik/N}\sum _{j=0}^{m-1}{y}_{2j+1}(e^{-2\pi ikj/m})\\ & =({\text{DFT}}_m({\mathbf{y}}_{\text{even}}){)}_k+{\omega }_N^{-k}\cdot ({\text{DFT}}_m({\mathbf{y}}_{\text{odd}}){)}_k\end{array}$$

Essentially, instead of solving one big problem of size $N$, we solve two half-sized problems (DFTs of size $N/2$), and then combine the results with only $O(N)$ extra steps (the “twiddle factors,” ${\omega }_N^{-k}$). By applying this recursive splitting many times ( ${\log }_{2}N$ times), the total cost comes down from $N^2$ to $N\log N$.

Applications and Practical Details

The FFT is not just a mathematical curiosity; it is a workhorse of modern science and engineering.

Frequency Analysis: The Signal’s DNA

The FFT is like a prism for data. It takes a complex signal and reveals the pure frequencies it’s composed of. The magnitude squared of the Fourier coefficients, $|c_k{|}^2$, is called the power spectrum. Peaks in the power spectrum correspond to the dominant frequencies in the signal.

Example: Finding Frequencies in a Noisy Signal

Let’s create a signal composed of two sine waves (at 1 Hz and 7 Hz) and add some random noise. Can we recover the original frequencies?

import numpy as np
import matplotlib.pyplot as plt
 
N = 128  # Number of sample points
t = np.linspace(0, 1, N, endpoint=False) # Time vector (samples taken over 1 second)
 
# Create a signal with two frequencies + noise
y = np.sin(1 * 2 * np.pi * t) + 0.5 * np.sin(7 * 2 * np.pi * t)
y += 0.5 * np.random.randn(N)
 
# Compute the DFT and the power spectrum
c = np.fft.fft(y)
power = np.abs(c)**2
 
# The frequencies corresponding to the coefficients (in Hz)
freqs = np.fft.fftfreq(N, d=1/N)
 
# Plotting
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.plot(t, y)
plt.title("Signal in Time Domain (Looks like noise!)")
plt.xlabel("Time")
 
plt.subplot(1, 2, 2)
# We only plot the first half of the frequencies, due to symmetry
plt.plot(freqs[:N//2], power[:N//2])
plt.title("Power Spectrum (The Signal's DNA)")
plt.xlabel("Frequency (Hz)")
plt.grid()
plt.show()

The plot of the power spectrum will show two sharp peaks at exactly 1 Hz and 7 Hz, rising clearly above the noise floor. The FFT has successfully extracted the underlying components from the seemingly random data.

Aliasing: The Wagon-Wheel Effect

A crucial question is: how many samples do we need? The Nyquist-Shannon sampling theorem provides the definitive answer. If the highest frequency present in a continuous signal is $f_{\max }$, you must sample at a rate greater than $2f_{\max }$ (the Nyquist rate) to avoid losing information.

If you sample too slowly, high frequencies get “folded down” and masquerade as lower frequencies. This is called aliasing. It’s the reason why in old movies, the wheels of a speeding stagecoach sometimes appear to spin slowly backward. The camera’s frame rate (the sampling frequency) is too slow to capture the rapid rotation of the wheel spokes, and our brain is aliased into seeing a slower, backward motion.

Beware of Aliasing

If you don’t sample fast enough, you will see ghosts. High frequencies will appear as low frequencies, and you will draw completely wrong conclusions from your data. Always be sure your sampling rate is at least twice the highest frequency you care about.

The Grand Finale: DFT and Chebyshev Interpolation

This chapter might seem far removed from the polynomials of Chapter 2, but here is the grand, unifying connection.

Recall the change of variables $x=\cos (\theta )$. This maps the interval $[-1,1]$ for $x$ (the standard polynomial domain) to the interval $[0,\pi ]$ for $\theta$.

• Chebyshev points in $x$ are of the form $\cos (k\pi /n)$. When you look at them in the $\theta$ variable, these are equally spaced points.
• A polynomial in $x$, $p(x)$, becomes the function $p(\cos \theta )$. A remarkable fact is that $p(\cos \theta )$ can be written as a sum of terms like $\cos (k\theta )$, which means it is a trigonometric polynomial in the $\theta$ variable.

This leads to a spectacular conclusion:

Interpolating a function $f(x)$ at Chebyshev nodes is mathematically identical to performing trigonometric interpolation on the related, periodic function $g(\theta )=f(\cos \theta )$ at equally spaced nodes in $\theta$.

This is a profound link between two worlds. It means we can use the incredible speed of the FFT to solve the polynomial interpolation problem. The problem of finding the coefficients for a Chebyshev expansion can be computed not in $O(n^2)$ time, but in $O(n\log n)$ time using a variant of the FFT called the Discrete Cosine Transform (DCT).

The stability and accuracy of Chebyshev polynomial interpolation are powered by the computational engine of the FFT. It’s a perfect example of how deep connections across different fields of mathematics lead to the most powerful and elegant numerical methods.

Next: Chapter 4 - Piecewise Polynomial Interpolation