Trigonometric interpolation and friends

The trigonometric interpolation problem

As usual, we start from a given interval \(I=[0, L)\) and \(N\) equidistant points \(x_k = k \frac{L}{N}\), \(k=0, 1, \ldots, N-1\). we now want to consider the following interpolation problem

Definition 15 (Trigonometric interpolation)

For a given list of sampled function values \(f_l = f(x_l)\), \(0\leqslant l \leqslant N-1\), find a (the?) trigonometric polynomial \(q(x)\) of the form

(46)\[ q_N(x) = \sum_{k=0}^{N-1} c_k e^{i 2 \pi k x/L} = \sum_{k=0}^{N-1} c_k \omega^{k}(x) \]

which interpolates the function values \(f_l\) at the points \(x_l\), i.e.

(47)\[ q_N(x_l) = f_l, \quad 0\leqslant l \leqslant N-1. \]

Similar to the Vandermonde matrix approach (sec:poly-interpol-direct) for the polynomial interpolation problem we discussed in Chapter Polynomial interpolation, we can set up a linear system for the coefficients \(c_k\),

\[ f_l = q_N(x_l) = \sum_{k=0}^{N-1} c_k \omega^{k}(x_l) = \sum_{k=0}^{N-1} c_k \omega_N^{kl}, \quad 0\leqslant l \leqslant N-1. \]

Or in matrix-vector notation

\[\begin{split} \begin{pmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega_N^{1} & \omega_N^{2} & \cdots & \omega_N^{(N-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega_N^{(N-1)} & \omega_N^{2(N-1)} & \cdots & \omega_N^{(N-1)^2} \end{pmatrix} \begin{pmatrix} c_0 \\ c_1 \\ \vdots \\ c_{N-1} \end{pmatrix} = \begin{pmatrix} f_0 \\ f_1 \\ \vdots \\ f_{N-1} \end{pmatrix}. \end{split}\]

Aha! This matrix is exactly the one we have seen in the definition on the inverse DFT, only that we now want to solve the inverse problem, i.e. we want to find the coefficients \(c_k\) for a given set of function values \(f_l\). But since the inverse DFT is invertible, and its inverse is given by the DFT matrix, we see that there is a unique solution for the coefficients \(c_k\). We record this observation in the following theorem.

Theorem 14 (Trigonometric interpolation problem)

There exists a unique trigonometric polynomial \(q(x)\) of the form (46) which solves the interpolation problem (47), and the coefficients \(c_k\) are given by the DFT of the sample vector \(\{f_l\}_{l=0}^{N-1}\), i.e.

\[ c_k = \widehat{f}_k = \frac{1}{N} \sum_{l=0}^{N-1} f_l \omega_N^{-l k}. \]

The trigonometric polynomial \(q(x)\) is called the trigonometric interpolant of the function values \(f_l\) at the points \(x_l\), and sometimes we write

\[ q_N(x) =: \pi^N f(x). \]

to emphasize the relation between the function \(f(x)\) and its trigonometric interpolant.

Often it desirable to compute a trigonmetric interpolating polynomial which frequencies are centered around \(0\).

More precisely, assuming that \(N = 2n + 1\) is odd, we wish to find a polynomial \(\tilde{q}_N(x)\) of the form

\[ \tilde{q}_N(x) = \sum_{k=-n}^{n} \tilde{c}_k \omega^{k}(x) \]

which satisfies the interpolation conditions \(\tilde{q}_N(x_l) = f_l\) for \(0\leqslant l \leqslant N-1\).

This one can in fact easily constructed from \(q_N(x)\) respectively the coefficients \(c_k\) as we will see now. We recall that the original polynomial \(q_N(x)\) satisfies

(48)\[\begin{align} f_l = q(x_l) &= \sum_{k=0}^{N-1} c_k \omega_N^{kl} \\ &= \sum_{k=0}^{n} c_k \omega_N^{kl} +\sum_{k=n+1}^{N-1} c_k \omega_N^{kl} \end{align}\]

for \( 0\leqslant l \leqslant N-1\).

Now we simply shift the index in the second sum by \(-N\) and use the periodicity of the \(N\)-th roots of unity \(\omega_N^{k+N} = \omega_N^k\) to see that

(49)\[\begin{align} q_N(x_l) &= \sum_{k=0}^{n} c_k \omega_N^{kl} +\sum_{k=-n}^{-1} c_{k+N} \omega_N^{(k+N)l} \\ &= \sum_{k=0}^{n} c_k \omega_N^{kl} +\sum_{k=-n}^{-1} c_{k+N} \omega_N^{k} = \sum_{k=-n}^{n} \widetilde{c}_k \omega_N^{kl} \end{align}\]

holds for \(0 \leqslant l \leqslant N-1\), where we set

\[\begin{split} \widetilde{c}_k = \begin{cases} c_k & \text{for } 0\leqslant k \leqslant n, \\ c_{k+N} & \text{for } -n \leqslant k \leqslant -1. \end{cases} \end{split}\]

Since \(\omega_N^{kl} = \omega^k(x_l)\) we thus conclude that the trigonometric polynomial

\[ \tilde{q}_N(x) = \sum_{k=-n}^{n} \widetilde{c}_k \omega^k(x) = \sum_{k=-(N-1)/2}^{(N-1)/2} \widetilde{c}_k \omega^k(x) \]

also satisfies the interpolation conditions \(\tilde{q}_N(x_l) = f_l\) for \(0\leqslant l \leqslant N-1\).

Note

This resembles more closely the form of the truncated complex Fourier series for an \(L\) periodic function \(f\),

\[ f \sim \sum_{k=-n}^{n} \widehat{f}(k) e^{i2\pi k x /L} = \sum_{k=-n}^{n} \widehat{f}(k) \omega^k(x) \]

For even \(N = 2n\), we can follow the same line of thoughts to see that with

\[\begin{split} \widetilde{c}_k = \begin{cases} c_k & \text{for } 0\leqslant k \leqslant n-1, \\ c_{k+N} & \text{for } -n \leqslant k \leqslant -1. \end{cases} \end{split}\]

the polynomial

\[ \tilde{q}(x) = \sum_{k=-n}^{n-1} \widetilde{c}_k \omega^k(x) \tilde{q}(x) = \sum_{k=-N/2}^{N/2-1} \widetilde{c}_k \omega^k(x) \]

satisfies the interpolation conditions \(\tilde{q}(x_l) = f_l\) for \(0\leqslant l \leqslant N-1\).

Because of the simple index shift, we can use the discrete fast Fourier transform of \(\{f_l\}_{l=0}^N\) to compute both \(c_k\) and \(\widetilde{c}_k\). As a matter of convention the following indexing is commonly used for the Fourier coefficients (But always check!):

For \(N = 2n +1 \) we have

(50)\[\begin{align} [\widehat{f}_0, \widehat{f}_1, \ldots, \widehat{f}_{N-1}] &= [c_0, c_1, \ldots, c_n, c_{n+1}, \ldots, c_{N-1}] \\ &= [ \widetilde{c}_0, \ldots, \widetilde{c}_n, \widetilde{c}_{-n}, \ldots, \widetilde{c}_{-1}] \\ &= [\widehat{f}_0, \widehat{f}_1, \ldots, \widehat{f}_{n}, \widehat{f}_{-n}, \ldots, \widehat{f}_{-1}] \end{align}\]

For \(N = 2n\) we have instead

(51)\[\begin{align} [\widehat{f}_0, \widehat{f}_1, \ldots, \widehat{f}_{N-1}] &= [c_0, c_1, \ldots, c_{n-1}, c_{n}, \ldots, c_{N-1}] \\ &= [ \widetilde{c}_0, \ldots, \widetilde{c}_{n-1}, \widetilde{c}_{-n}, \ldots, \widetilde{c}_{-1}] \\ &= [\widehat{f}_0, \widehat{f}_1, \ldots, \widehat{f}_{n-1}, \widehat{f}_{-n}, \ldots, \widehat{f}_{-1}] \end{align}\]

Real trigonometric interpolation

Similar to the different representations of the Fourier series, cf. fou:eq:fourier-partial-real, we sometimes want to represent the trigonometric interpolant in terms of \(\cos\) and \(\sin\) functions. This is particularly useful to arrive at real-valued trigonometric interpolant in the case that the sampled function real-valued.

To arrive at such a rewrite representation, we proceed exactly in the same way as in the case of the Fourier series: Simply rewrite the complex exponential functions in terms of \(\cos\) and \(\sin\) functions via Euler’s formula and collect the terms with the same trigonometric functions.

Let’s start with the case of an odd \(N = 2n + 1\). Then

(52)\[\begin{align} \widetilde{q}_N(x) = \sum_{k=-n}^{n} c_k e^{i 2 \pi k x/L} &= c_0 + \sum_{k=1}^{n} c_k e^{i 2 \pi k x/L} + c_{-k} e^{-i 2 \pi k x/L} \\ &= c_0 + \sum_{k=1}^{n} c_k \left( \cos(2 \pi k x/L) + i \sin(2 \pi k x/L) \right) \\ &\phantom{= c_0} + \sum_{k=1}^{n} c_{-k} \left( \cos(2 \pi k x/L) - i \sin(2 \pi k x/L) \right) \\ &= c_0 + \sum_{k=1}^{n} \left( (c_k + c_{-k}) \cos(2 \pi k x/L) + i(c_k - c_{-k}) \sin(2 \pi k x/L) \right) \\ &= \dfrac{a_0}{2} + \sum_{k=1}^{n} \left( a_k \cos(2 \pi k x/L) + b_k \sin(2 \pi k x/L) \right) \end{align}\]

where we set \(a_k = c_k + c_{-k}\) and \(b_k = i(c_k - c_{-k})\).

For real-valued samples \(\{f_l\}_{l=0}^{N-1}\), we can now make the following observation:

Since \(\overline{f_l} = f_l\) for real-valued \(f_l\), we infer that \(\overline{c_k} = c_{-k}\) since

\[ \overline{c_k} = \overline{\widehat{f}_k} = \overline{\frac{1}{N} \sum_{l=0}^{N-1} f_l \omega_N^{-lk}} = \dfrac{1}{N} \sum_{l=0}^{N-1} f_l \omega_N^{lk} = \widehat{f}_{-k} = c_{-k} \]

Consequently,

\begin{align} a_k &= c_k + c_{-k} = c_k + \overline{c_k} = 2 \Re(c_k) = 2 \Re(\widehat{f}_k) = \frac{2}{N} \Re \left( \sum_{l=0}^{N-1} f_l \omega_N^{-lk} \right) \\ &= \frac{2}{N} \sum_{l=0}^{N-1} f_l \cos(2 \pi k l/N) \frac{2}{N} = \sum_{l=0}^{N-1} f_l \cos(2 \pi k x_l/L) \end{align}

Similarly, we have that

(54)\[\begin{align} b_k &= i(c_k - c_{-k}) = i(c_k - \overline{c_k}) = 2 \Im(c_k) = 2 \Im(\widehat{f}_k) = \frac{2}{N} \Im \left( \sum_{l=0}^{N-1} f_l \omega_N^{-lk} \right) \\ &= \frac{2}{N} \sum_{l=0}^{N-1} f_l \sin(2 \pi k l/N) = \frac{2}{N} \sum_{l=0}^{N-1} f_l \sin(2 \pi k x_l/L) \end{align}\]

Consequently, \(a_k\) and \(b_k\) are real-valued, and the trigonometric interpolant \(\widetilde{q}_N(x)\) is in fact a real-valued trigonometric polynomial for all \(x\in [0,L)\) (and not only for the sampled points \(x_l\))!

For even \(N = 2n\), we follow the convention above and we consider again

(55)\[\begin{align} \widetilde{q}_N(x) = \sum_{k=-n}^{n-1} c_k e^{i 2 \pi k x/L} = c_{-n} e^{-i 2 \pi n x/L} + \sum_{k=-(n-1)}^{n-1} c_k e^{i 2 \pi k x/L} \end{align}\]

We can rewrite the sum in the second term as above to obtain

\[ \widetilde{q}_N(x) = c_{-n} e^{- i 2 \pi n x/L} + \dfrac{a_0}{2} + \sum_{k=1}^{n-1} \left( a_k \cos(2 \pi k x/L) + b_k \sin(2 \pi k x/L) \right) \]

As before, we conclude that \(a_k\) and \(b_k\) are real-valued if \(\{f_l\}_{l=0}^{N-1}\) are real-valued, and so is the resulting superposition of \(\cos\) and \(\sin\) functions.

But what about the term \(c_{-n} e^{- i 2 \pi n x/L}\)?

First let us observe that the term \(c_{-n}\) is in fact real-valued, since (recalling that \(N = 2n\))

\[ c_{-n} = \dfrac{1}{N} \sum_{l=0}^{N-1} f_l e^{i 2 \pi n l/(2n)} = \dfrac{1}{N} \sum_{l=0}^{N-1} f_l (-1)^l = \dfrac{1}{N} \sum_{l=0}^{N-1} f_l \cos(2\pi n l/N)= \dfrac{1}{N} \sum_{l=0}^{N-1} f_l \cos(2\pi n x_l/L) \]

In particular this means, that \(\widetilde{q}_N(x)\) is not a purely real-valued trigonometric polynomial, it only happens to be real-valued at the sample points \(x_l\)!

To remedy this situation and obtain a real-valued trigonometric interpolant, we recall that

\[ \Re(\widetilde{q}_N(x_l)) = \widetilde{q}_N(x_l) \]

must hold for all \(0\leqslant l \leqslant N-1\) since \(\widetilde{q}_N(x_l) = f_l\) is real-valued. Consequently, we can simply \(\Re(\widetilde{q}_N(x))\) to find a real-valued trigonometric interpolating polynomial.

But since \(c_{-n}\) is real-valued, and thus \(\Re(c_{-n} e^{- i 2 \pi n x/L}) = c_{-n} \Re( e^{- i 2 \pi n x/L})\), we see that (after setting \(2 a_n := c_{-n}\))

(56)\[\begin{align} \Re(\widetilde{q}_N(x)) &= \Re(\tfrac{a_n}{2} e^{- i 2 \pi n x/L}) + \dfrac{a_0}{2} + \sum_{k=1}^{n-1} \left( a_k \cos(2 \pi k x/L) + b_k \sin(2 \pi k x/L) \right) \\ &= \dfrac{a_n}{2} \cos(2 \pi n x/L) + \dfrac{a_0}{2} + \sum_{k=1}^{n-1} \left( a_k \cos(2 \pi k x/L) + b_k \sin(2 \pi k x/L) \right) \end{align}\]

is a real-valued trigonometric interpolant of the real-valued function values \(f_l\) at the points \(x_l\).

Let us summarize this in the following theorem.

Theorem 15 (Real trigonometric interpolation)

Let \(\{f_l\}_{l=0}^{N-1}\) be a list of real-valued function values at the points \(x_l = l \frac{L}{N}\), \(0\leqslant l \leqslant N-1\). Define the coefficients \(a_k\) and \(b_k\) by

\[ a_k = \sum_{l=0}^{N-1} f_l \cos(2 \pi k x_l/L), \quad b_k = \sum_{l=0}^{N-1} f_l \sin(2 \pi k x_l/L). \]

If \(N = 2n + 1\) is odd, then the real-valued trigonometric polynomial \(p_n(x)\) of order \(n\) given by

(57)\[ p_n(x) = \dfrac{a_0}{2} + \sum_{k=1}^{n} \left( a_k \cos(2 \pi k x/L) + b_k \sin(2 \pi k x/L) \right) \]

satisfies \(f_l = p_n(x_l)\) for \(0\leqslant l \leqslant N-1\).

If \(N = 2n\) is even on the other hand, then the real-valued interpolating trigonometric polynomial \(p_n(x)\) of order \(n\) given by

(58)\[ p_n(x) = \dfrac{a_0}{2} + \sum_{k=1}^{n-1} \left( a_k \cos(2 \pi k x/L) + b_k \sin(2 \pi k x/L) \right) + \dfrac{a_n}{2} \cos(2 \pi n x/L). \]

Best approximation properties

We now take another at the properties of the discrete Fourier transform, the resulting coefficients \(\widehat{f}_k\) and the properties of trigonometric polynomials constructed from them.

As usual we assume we have a function \(f(x)\) which is \(L\)-periodic and we sample it at \(N\) equidistant points \(x_l = l \frac{L}{N}\), \(0\leqslant l \leqslant N-1\).

We have seen that the trigonometric interpolant

\[ \pi^N f(x) := q_N(x) := \sum_{k=0}^{N-1} \widehat{f}_k \omega^k(x) \]

interpolates the function values \(f_l = f(x_l)\) at the points \(x_l\).

But what happens if we truncate this sum and consider a lower order trigonometric polynomial of the form

\[ \pi^N_m f(x) := q_m(x) := \sum_{k=0}^{m} \widehat{f}_k \omega^k(x) \]

for \(0\leqslant m < N-1\)?

To gain some insight into this question, we change slightly our perspective and focus on the vector presentation of the sampled function values \(\{f_l\}_{l=0}^{N-1}\). We write

\[ \mathbf{F} = [f_0, f_1, \ldots, f_{N-1}] \]

As we noted before, the \(N\) vectors

\[ \mathbf{W^k} = (\omega^k(x_l))_{l=0}^{N-1} = (\omega_N^{lk})_{l=0}^{N-1} = [1, \omega_N^k, \omega_N^{2k}, \ldots, \omega_N^{(N-1)k}] \quad {k=0, 1, \ldots, N-1} \]

define a orthonormal basis of \(\mathbb{C}^N\) with respect to the discrete inner product \(\langle \cdot, \cdot \rangle_N\).

Now the interpolation property \(\pi^N f(x_l) = q_N(x_l) = f_l\) for \(0\leqslant l \leqslant N-1\) is nothing else than the statement that the vector \(\mathbf{F}\) can be written as a linear combination of the vectors \(\mathbf{W^l}\), and the coefficients of this linear combination are given by the DFT of \(\mathbf{F}\).

(59)\[\begin{align} \mathbf{F} &= [f_0, f_1, \ldots, f_{N-1}] \\ &= [q_N(x_0), q_N(x_1) \ldots, q_N(x_{N-1})] \\ &= \left[ \sum_{k=0}^{N-1} \widehat{f}_k \omega^k(x_0), \ldots, \sum_{k=0}^{N-1} \widehat{f}_k \omega^k(x_{N-1}) \right] \\ &= \sum_{k=0}^{N-1} \widehat{f}_k \left[ \omega^k(x_0), \ldots, \omega^k(x_{N-1}) \right] \\ &= \sum_{k=0}^{N-1} \widehat{f}_k \left[ 1, \omega_N^k, \ldots, \omega_N^{k(N-1)} \right] \\ &= \sum_{k=0}^{N-1} \widehat{f}_k \mathbf{W^k} \end{align}\]

But do you remember this definition of the \(\widehat{f}_k\)? Recall that we could write them with the discrete inner product as

(60)\[\begin{align} \widehat{f}_k &= \langle f(x), \omega^k(x) \rangle_N \\ &= \langle (f_l)_{l=0}^{N-1}, (\omega^l(x_l))_{l=0}^{N-1} \rangle_N \\ &= \langle \mathbf{F}, \mathbf{W^k} \rangle_N \end{align}\]

In other words

\[ \mathbf{F} = \sum_{k=0}^{N-1} \langle \mathbf{F}, \mathbf{W}^k \rangle_N \mathbf{W^k} \]

Note that this exactly the usual way for computing the presentation of a vector in terms of an orthonormal basis.

It we repeat this argument for the truncated trigonometric polynomial \(\pi^N_m f(x) = q_m(x) \sum_{k=0}^{m} \widehat{f}_k \omega^k(x)\), we obtain

(61)\[\begin{align} \mathbf{Q}_m &:= [q_m(x_0), q_m(x_1) \ldots, q_m(x_{N-1})] \\ &= \left[ \sum_{k=0}^{m} \widehat{f}_k \omega^k(x_0), \ldots, \sum_{k=0}^{m} \widehat{f}_k \omega^k(x_{N-1}) \right] \\ &= \sum_{k=0}^{m} \widehat{f}_k \mathbf{W^k} \\ &= \sum_{k=0}^{m} \langle \mathbf{F}, \mathbf{W}^k \rangle_N \mathbf{W^k} \end{align}\]

Aha! We see that \(\mathbf{Q}_m\) is nothing else than the orthogonal projection \(\Pi_m \mathbf{F}\) of \(\mathbf{F}\) onto the \(m\) dimensional subspace of \(\mathbb{C}^N\),

\[ V_m = \mathrm{span}(\{\mathbf{W^0}, \ldots, \mathbf{W}^m\}) \]

with respect to the discrete inner product \(\langle \cdot, \cdot \rangle_N\)!

Equivalently, we can say that the truncated trigonometric polynomial

(62)\[\begin{align} \pi^N_m f(x) = q_m(x) &= \sum_{k=0}^{m} \widehat{f}_k \omega^k(x) \\ &= \sum_{k=0}^{m} \langle f, \omega^k\rangle_N \omega^k(x) \end{align}\]

is the orthogonal projection of the function \(f(x)\) onto the subspace of trigonometric polynomials of order \(m\)

\[ T_m^{\mathbb{C}} = \mathrm{span}(\{\omega^0, \ldots, \omega^m\}) = \{ r_m(x) = \sum_{k=0}^{m} d_k \omega^k(x) \mid d_k \in \mathbb{C}\} \]

with respect to the discrete inner product \(\langle \cdot, \cdot \rangle_N\).

Now remember that in general, for orthogonal projection \(\Pi_m \mathbf{F}\) of a vector \(\mathbf{F}\) onto a subspace \(V_m\), the projection error \(\mathbf{F}- \Pi_m \mathbf{F}\) is orthogonal to the subspace \(V_m\),

(63)\[ \langle \mathbf{F}- \Pi_m \mathbf{F}, \mathbf{R}_m \rangle_N = 0 \quad \text{for all } \mathbf{R}_m \in V_m, \]

or equivalently

(64)\[ \langle f - \pi^N_m f, r_m \rangle_N = 0 \quad \text{for } r_m \in T_m^{\mathbb{C}}. \]

We can use this property to derive a best approximation property of the truncated trigonometric polynomial \(\pi^N_m f(x)\).

Theorem 16 (Best approximation property of the truncated trigonometric polynomial)

The truncated trigonometric polynomial \(\pi^N_m f(x)\) satisfies a best approximation property in the sense that

\[ \| f - \pi^N_m f \|_N = \min_{r_m \in T_{\mathbb{C}}^N} \| f - r_m \|_N. \]

In other words, \(\pi^N_m f(x)\) minimizes the squared error

\[ \sum_{l=0}^{N-1} |f_l - r_m(x_l)|^2 \]

among all trigonometric polynomials \(r_m\) of order \(m\). We say it say \(\pi^N_m f(x)\) has the least square error. among all trigonometric polynomials of order \(m\).

Proof.

For any given trigonometric polynomial \(r_m(x) = \sum_{k=0}^{m} d_k \omega^k(x)\) of order \(m\), we use the orthogonality property (64) to see that

(65)\[\begin{align} \| f - \pi^N_m f \|_N^2 &= \langle f - \pi^N_m f, f - \pi^N_m f \rangle_N \\ &= \langle f - \pi^N_m f, f - r_m + r_m - \pi^N_m f \rangle_N \\ &= \langle f - \pi^N_m f, f - r_m \rangle_N + \underbrace{\langle f - \pi^N_m f, r_m - \pi^N_m f \rangle_N}_{=0 \text{ by orthogonality}} \\ &= \langle f - \pi^N_m f, f - r_m \rangle_N \\ & \leqslant \| f -\pi^N_m f \|_N \| f - r_m \|_N \end{align}\]

Assuming If \(\| f -\pi^N_m f \|_N = 0\)

Otherwise if \(\| f -\pi^N_m f \|_N > 0\), we can divide by this term to obtain

\[ \| f - \pi^N_m f \|_N \leqslant \| f - r_m \|_N \]

Since \(r_m\) was an arbitrary trigonometric polynomial of order \(m\), we see that the truncated trigonometric polynomial \(\pi^N_m f(x)\) satisfies the minimization property

\[ \| f - \pi^N_m f \|_N = \min_{r_m \in T_{\mathbb{C}}^N} \| f - r_m \|_N. \]