Fourier transform - Numerical Analysis

#%config InlineBackend.figure_format = 'svg'
from pylab import *

Periodic functions can be approximated by a series of sine and cosine functions. This leads us to the Fourier transform of a function defined on a bounded interval and we can think of this function being extended periodically to $\re$ .

6.11.1Fourier series¶

Let $f : [-\pi,\pi] \to \re$ be a function which is extended periodically to the whole real line. The Fourier series is defined as

Sf(x) = a_0 + \sum_{k=1}^\infty [ a_k \cos(k x) + b_k \sin(k x) ]

(1)

where

\begin{align*} a_0 &= \frac{1}{2\pi} \int_{-\pi}^\pi f(x) \ud x \\ a_k &= \frac{1}{\pi} \int_{-\pi}^\pi f(x) \cos(k x) \ud x, \qquad k=1, 2, \ldots \\ b_k &= \frac{1}{\pi} \int_{-\pi}^\pi f(x) \sin(k x) \ud x, \qquad k=1, 2, \ldots \end{align*}

(2)

Due to periodicity $k$ only takes integer values.

The Fourier series can be written as

Sf(x) = \sum_{k=-\infty}^\infty \hat f_k \ee^{\ii k x}, \qquad \hat f_k = \frac{1}{2\pi} \int_{-\pi}^\pi f(x) \ee^{-\ii k x} \ud x

(3)

where

\hat f_0 = a_0, \qquad \hat f_k = \begin{cases} \half(a_k + \ii b_k), & k = -1, -2, \ldots \\ \half(a_k - \ii b_k), & k = +1, +2, \ldots \end{cases}

(4)

If $f(x)$ is real, $a_k, b_k$ are real and we have

\hat f_{-k} = \conj{\hat f_k}, \qquad k=1,2,\ldots

(5)

6.11.2Convergence of Fourier series¶

When is $f = Sf$ ? To talk of such questions, let us define the truncated series

S_n f(x) = a_0 + \sum_{k=1}^n [ a_k \cos(k x) + b_k \sin(k x) ] = \sum_{k=-n}^n \hat f_k \ee^{\ii k x}

(8)

Note that we have an expansion in terms of an orthogonal basis $\phi_k(x) = \ee^{\ii k x}$ , since

\int_{-\pi}^\pi \phi_j(x) \phi_k(x) \ud x = 2 \pi \delta_{jk}

(9)

There are three types of convergence wrt $n \to \infty$ .

pointwise convergence
$S_n f(x) \to f(x), \qquad \forall x$
(10)
uniform convergence
$\norm{S_n f - f}_\infty \to 0$
(11)
mean square convergence
$\norm{S_n f - f}_2 \to 0$
(12)

where the norms are defined as

\norm{f}_\infty = \max_{x \in [-\pi,\pi]} |f(x)|, \qquad \norm{f}_2 = \left( \int_{-\pi}^\pi |f(x)|^2 \ud x \right)^\half

(13)

Uniform convergence is the strongest property, and it implies pointwise and mean square convergence Tveito & Winther, 2005.

Clearly we need $|\hat f_k| \to 0$ as $|k| \to \infty$ for the series to converge and the rate of decay needs to be sufficiently fast. The rate of decay of $\hat f_k$ depends on the smoothness of the function.

Definition 1

The function $f : [-\pi,\pi] \to \re$ is said to be piecewise continuous if it is continuous at all but a finite number of points where both left and right limits exist.
The function $f : [-\pi, \pi] \to \re$ is said to be piecewise differentiable if it is differentiable everywhere except for a finite number of points where the one-sided derivatives exist.
We say that $f \in \cts_p^0[-\pi,\pi]$ if it is continuous and $f(-\pi) = f(\pi)$ ; i.e., the periodic extension of $f$ to $\re$ is continuous.
We say that $f \in \cts_p^\nu[-\pi,\pi]$ if it is ν-times continuously differentiable and
$f^{(s)}(-\pi) = f^{(s)}(\pi), \qquad s = 0,1,\ldots,\nu$
(14)
i.e., the periodic extension of $f$ to $\re$ is in $\cts^\nu(\re)$ .

The term piecewise continuous can be replaced with bounded variation. The proof is based on an adaptation of Gander & Kwok, 2018, Theorem 4.2.

Proof 1

(1a) First consider $f(x)$ to be a step function

f(x) = y_j, \qquad x \in (x_{j-1},x_j), \qquad 1 \le j \le m

(17)

for some partition $\{x_0, x_1, \ldots, x_m\}$ of $[-\pi,\pi]$ . Then

\begin{align} \hat f_k &= \frac{1}{2\pi} \int_{-\pi}^\pi f(x) \ee^{-\ii k x} \ud x \\ &= \frac{1}{2\pi \ii k} \left[ y_1 + \ee^{-\ii k x_1} (y_2 - y_1) + \ee^{-\ii k x_2} (y_3 - y_2) + \right. \\ & \qquad\qquad \left. \ldots + \ee^{-\ii k x_{n-1}} (y_m - y_{m-1}) - y_m \right] \end{align}

(18)

and hence

|\hat f_k| \le \frac{\TV(f) + |f(\pi) - f(-\pi)|}{2 \pi |k|} = \order{ \frac{1}{|k|} }

(19)

where

\TV(f) = \sum_{j=2}^m |y_j - y_{j-1}| < \infty

(20)

(1b) Now let $f(x)$ be a general piecewise continuous function; it is Riemann integrable and for any $\epsilon > 0$ , there is a partition $P = \{ x_0, x_1, \ldots, x_m \}$ such that (Rudin, 1976, Theorem 6.6)

U(f,P) - L(f,P) = \sum_{j=1}^m M_j (x_j - x_{j-1}) - \sum_{j=1}^m m_j (x_j - x_{j-1}) \le \epsilon

(21)

where

m_j = \min_{x \in [x_{j-1},x_j]} f(x) = f(\xi_j), \qquad M_j = \max_{x \in [x_{j-1},x_j]} f(x)

(22)

Consider the step function

\tilde f(x) = m_j, \qquad x \in (x_{j-1}, x_j)

(23)

for which

\TV(\tilde f) = \sum_{j=2}^m |f(\xi_j) - f(\xi_{j-1})| \le \TV(f)

(24)

Now (Rudin, 1976, Theorem 6.12, 6.13)

\begin{align} \left| \int_{-\pi}^\pi f(x) \ee^{-\ii k x} \ud x - \int_{-\pi}^\pi \tilde f(x) \ee^{-\ii k x} \ud x \right| &\le \int_{-\pi}^\pi |f(x) - \tilde f(x)| \ud x \\ &= \sum_{j=1}^m \int_{x_{j-1}}^{x_j} |f(x) - \tilde f(x)| \ud x \\ &= \sum_{j=1}^m \int_{x_{j-1}}^{x_j} (f(x) - m_j) \ud x \\ &\le \sum_{j=1}^m \int_{x_{j-1}}^{x_j} (M_j - m_j) \ud x \\ &= U(f,P) - L(f,P) \\ &\le \epsilon \end{align}

(25)

Then the Fourier transform

\begin{align} \hat f_k &= \frac{1}{2\pi} \int_{-\pi}^\pi f(x) \ee^{-\ii k x} \ud x \\ &\approx \frac{1}{2\pi} \int_{-\pi}^\pi \tilde f(x) \ee^{-\ii k x} \ud x \\ &\le \frac{\TV(\tilde f) + |\tilde f(\pi) - \tilde f(-\pi)|}{2 \pi |k|} \qquad \textrm{using Part (1a)} \\ &\lessapprox \frac{\TV(f) + |f(\pi) - f(-\pi)|}{2 \pi |k|} \\ &= \order{\frac{1}{|k|}} \end{align}

(26)

(2) If $f \in \cts^0_p[-\pi,\pi]$ and piecewise differentiable, i.e., $s=1$ , then using integration by parts

\begin{align} \hat f_k &= \frac{1}{2\pi} \int_{-\pi}^\pi f(x) \ee^{-\ii k x} \ud x \\ &= -\frac{1}{2 \pi \ii k} \left[ f(x) \ee^{-\ii k x} \right]_{-\pi}^\pi + \frac{1}{\ii k} \left[ \frac{1}{2 \pi} \int_{-\pi}^\pi f'(x) \ee^{-\ii k x} \ud x \right] \\ & \qquad \textrm{Using the result of Part (1) in the second term above} \\ &= -\frac{\cos(k\pi)}{2 \pi \ii k} [f(\pi) - f(-\pi)] + \frac{1}{\ii k} \order{\frac{1}{k}} \\ &= \order{\frac{1}{k^2}}, \qquad \textrm{since } f(-\pi) = f(\pi) \end{align}

(27)

We see that every time we do an integration by parts, we get a factor of $1/k$ . For any $s > 1$ , we can perform several integration by parts and reach the conclusion of Part (2).

Example 1 (Piecewise continuous function)

f(x) = \begin{cases} 1, & \half\pi < x \le \frac{3}{2}\pi \\ 0, & \textrm{otherwise} \end{cases}, \qquad x \in [0,2\pi]

(28)

f = lambda x: (x > 0.5*pi)*(x < 1.5*pi)*(1)
x = linspace(0,2*pi,500)
plot(x,f(x))
grid(True), xlabel('x'), ylabel('f(x)');

Its Fourier coefficients are

\hat f_k = \begin{cases} \half, & k = 0 \\ 0, & k \textrm{ even} \\ \frac{(-1)^{(k+1)/2}}{\pi k}, & k \textrm{ odd} \end{cases}, \qquad |\hat f_k| = \order{\frac{1}{|k|}}

(29)

This corresponds to Part (1) of Theorem 1.

Example 2 (Continuous function)

f(x) = \begin{cases} 2x/\pi - 1, & \half\pi \le x \le \pi \\ 3 - 2x/\pi, & \pi \le x \le \frac{3}{2}\pi \\ 0, & \textrm{otherwise} \end{cases}, \qquad x \in [0,2\pi]

(30)

f = lambda x: (x > 0.5*pi)*(x < pi)*(2*x/pi - 1) \
              + (x > pi)*(x < 1.5*pi)*(3-2*x/pi)
x = linspace(0,2*pi,500)
plot(x,f(x))
grid(True), xlabel('x'), ylabel('f(x)');

|\hat f_k| = \order{ \frac{1}{|k|^2} }

(31)

This corresponds to Part (2) of Theorem 1 with $s=1$ .

Example 3 (Infinitely differentiable but not periodic)

f(x) = \sin(x/2), \qquad x \in [0,2\pi]

(32)

f = lambda x: sin(x/2)
df= lambda x: cos(x/2)/2
x = linspace(0,2*pi,500)
plot(x,f(x),label='f(x)')
plot(x,df(x),label='f\'(x)')
legend(), grid(True), xlabel('x');

Its Fourier coefficients are

\hat f_k = \frac{2}{\pi(1 - 4 k^2)}, \qquad |\hat f_k| = \order{\frac{1}{k^2}}, \qquad |k| \to \infty

(33)

Now $f(0) = f(2\pi)$ but

f'(0) = \half \cos(0) = \half \qquad\ne\qquad f'(2\pi) = \half \cos(\pi) = -\half

(34)

and hence $f \in \cts^0_p[0,2\pi]$ but $f \notin \cts^1_p[0,2\pi]$ . This corresponds to Part (1) of Theorem 1 with $s=1$ .

We can see this with integration-by-parts;

\begin{align} \hat f_k &= \frac{1}{2\pi} \int_{0}^{2\pi} f(x) \ee^{-\ii k x} \ud x \\ &= \frac{1}{2 \pi \ii k} \int_{0}^{2\pi} f'(x) \ee^{-\ii k x} \ud x \qquad \textrm{since } f(0) = f(2\pi) \\ &= \order{\frac{1}{k^2}} [ f'(2\pi) - f'(0)] + \frac{1}{2 \pi (\ii k)^2} \int_{0}^{2\pi} f''(x) \ee^{-\ii k x} \ud x \end{align}

(35)

Though we can perform more integrations giving rise to terms of order $1/|k|^3, 1/|k|^4, \ldots$ ,

|\hat f_k| = \order{\frac{1}{k^2}} [ f'(2\pi) - f'(0)] + \order{\frac{1}{|k|^3}} + \order{\frac{1}{|k|^4}} + \ldots

(36)

the size of $|\hat f_k|$ is $\order{1/k^2}$ for $|k| \to \infty$ , since $f'(0) \ne f'(2\pi)$ .

Remark 2

We can use Sympy to compute the Fourier transform. This code is for the above function.

def compute_fk():
    from sympy import pi, sin, exp, Symbol, integrate, pprint, simplify
    k = Symbol('k', integer=True)
    x = Symbol('x', real=True)
    fk = 1/(2*pi) * integrate(sin(x/2)*exp(-1j*k*x),(x,0,2*pi))
    pprint(simplify(fk))

compute_fk()

      -0.5       
─────────────────
  ⎛     2       ⎞
π⋅⎝1.0⋅k  - 0.25⎠

Example 4 (Infinitely differentiable and periodic)

f(x) = \frac{3}{5 - 4\cos x}, \qquad x \in [0,2\pi]

(37)

f = lambda x: 3.0/(5.0 - 4.0*cos(x))
x = linspace(0,2*pi,500)
plot(x,f(x))
grid(True), xlabel('x'), ylabel('f(x)');

Its Fourier coefficients are

\hat f_k = \frac{1}{2^{|k|}}

(38)

$f \in \cts^s_p[0,2\pi]$ for every $s \ge 0$ , and the Fourier coefficients decay faster than any power of $1/|k|$ , which is consistent with exponential decay.

Theorem 2 (Convergence of Fourier series)

(1) If $f : [-\pi,\pi] \to \re$ is piecewise continuous, then

\lim_{n \to \infty} \norm{S_n f - f}_2 = 0

(39)

Moreover, Parseval’s relation holds

\sum_{k=-\infty}^\infty |\hat f_k|^2 = \norm{f}_2^2

(40)

(2) If $f \in \cts_p^0[-\pi,\pi]$ and $f'$ is piecewise continuous, then the Fourier series converges absolutely and uniformly,

\lim_{n \to \infty} \norm{S_n f - f}_\infty = 0

(41)

Moreover, the Fourier transform of $g=f'$ is

\hat g_k = \ii k \hat f_k

(42)

Example 5 (Continuation of Example 1)

Here are the truncated Fourier series for the square hat function.

def fn(x,n):
    f = 0.5 * ones_like(x)
    for k in range(1,n,2):
        e1 = (-1)**(( k+1)//2) * exp( 1j*k*x) / ( k*pi)
        e2 = (-1)**((-k+1)//2) * exp(-1j*k*x) / (-k*pi)
        f += real(e1 + e2)
    return f

x = linspace(0,2*pi,1000)
f = (x > 0.5*pi)*(x < 1.5*pi)*(1) + 0.0
for i,n in enumerate([10,20,30,50]):
    subplot(2,2,i+1)
    plot(x, f, x, fn(x,n))
    text(pi,0.5,'n = ' + str(n),ha='center')
    grid(True);

As $n$ increases, we see that the errors do not decrease around the discontinuities and there is no convergence in maximum norm, but there is convergence in 2-norm as stated in Part (1) of Theorem 2.

We have observed this kind of Gibbs oscillations when we interpolate a discontinuous function with polynomials.

Moreover, it looks like there is pointwise convergence away from the discontinuities, and this is the case, see next Theorem.

Exercise 2

Compute the Fourier transform of the triangular hat function in Example 2. Verify the Fourier transform of the functions in the other examples also.

Plot the truncated Fourier series $(S_n f)(x)$ of all the functions in the above examples for a few values of $n$ as we did for the square hat function.

To test convergence as $n \to \infty$ , compute $S_n f$ on a uniform grid $\{ x_j \}$ of say $m = 1000$ points and measure the error as

\begin{align} \norm{S_n f - f}_\infty &\approx \max_{j} |S_n f(x_j) - f(x_j)| \\ \norm{S_n f - f}_2 &\approx \left( \frac{2\pi}{m} \sum_{j} [S_n f(x_j) - f(x_j)]^2 \right)^\half \end{align}

(44)

Plot error versus $n$ in loglog scale or semilogy scale; choose the scale according to the smoothness of the function.

References¶

Tveito, A., & Winther, R. (2005). Introduction to Partial Differential Equations (Vol. 29). Springer-Verlag. 10.1007/b138016
Davis, P. J. (1963). Interpolation and Approximation. Dover Publications.
Shen, J., Tang, T., & Wang, L.-L. (2011). Spectral Methods: Algorithms, Analysis and Applications (Vol. 41). Springer Berlin Heidelberg. 10.1007/978-3-540-71041-7
Gander, M. J., & Kwok, F. (2018). Numerical Analysis of Partial Differential Equations Using Maple and MATLAB. Society for Industrial. 10.1137/1.9781611975314
Rudin, W. (1976). Principles of Mathematical Analysis (3rd ed.). McGraw Hill.

Numerical Analysis

Spline functions

Numerical Analysis

Trigonometric interpolation

6.11Fourier transform

6.11.1Fourier series¶

6.11.2Convergence of Fourier series¶