BVH:MC:Author needs to add figure labels

Import standard modules:


In [ ]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from IPython.display import HTML 
HTML('../style/course.css') #apply general CSS

Import section specific modules:


In [ ]:
from IPython.display import HTML
from ipywidgets import interact
HTML('../style/code_toggle.html')

2.3. Fourier Series

While Fourier series are not immediately required to understand the required calculus for this book, they are closely connected to the Fourier transform, which is an essential tool. Moreover, we noticed a few times that the principle of the harmonic analysis or harmonic decomposition is essential- despite its simplicity - often not fully understood. We hence give a very brief summary, not caring about existence questions.

2.3.1 Definition

The Fourier series of a function $f: \mathbb{R} \rightarrow \mathbb{R}$ with real coefficients is defined as

$$ f_{\rm F}(x) \,=\, \frac{1}{2}c_0+\sum_{m = 1}^{\infty}c_m \,\cos(mx)+\sum_{m = 1}^{\infty}s_m \,\sin(mx), $$

with the Fourier coefficients $c_m$ and $s_m$

$$ \left( c_0 \,=\,\frac{1}{\pi}\int_{-\pi}^{\pi}f(x)\,dx \right)\\ c_m \,=\,\frac{1}{\pi}\int_{-\pi}^{\pi}f(x)\,\cos(mx)\,dx \qquad m \in \mathbb{N_0}\\ s_m \,=\,\frac{1}{\pi}\int_{-\pi}^{\pi}f(x)\,\sin(mx)\,dx \qquad m \in \mathbb{N_0}.\\{\rm } $$

If $f_{\rm F}$ exists, it is identical to $f$ in all points of continuity. For functions which are periodic with a period of $2\pi$ the Fourier series converges. Hence, for continuous periodic function with a period of $2\pi$ the Fourier series converges and $f_{\rm F}=f$.

The Fourier series of a function $f: \mathbb{R} \rightarrow \mathbb{R}$ with imaginary coefficients is defined as

$$ f_{\rm IF}(x) \,=\, \sum_{m = -\infty}^{\infty}a_m \,e^{\imath mx},$$

with the Fourier coefficients $a_m$

$$ a_m \,=\, \frac{1}{2\pi}\int_{-\pi}^{\pi}f(x)e^{-\imath mx}\,dx\qquad \forall m\in\mathbb{Z}.$$

The same convergence criteria apply and one realisation can be transformed to the other. Making use of Euler's formula ➞ , one gets

$$ \begin{split} a_m \,&=\, \frac{1}{2\pi}\int_{-\pi}^{\pi}f(x)\,[\cos(mx)-\imath \,\sin(mx)]\,dx &=\,\left \{ \begin{array}{lll} \frac{1}{2} (c_m+\imath s_m) & {\rm for} & m < 0\\ \frac{1}{2} c_m & {\rm for} & m = 0\\ \frac{1}{2} (c_m-\imath\,s_m) & {\rm for} & m > 0\\ \end{array} \right. \end{split}, $$

and accordingly, $\forall m \in \mathbb{N_0}$,

$$ \\ \begin{split} c_m \,&=\, a_m+a_{-m}\\ s_m \,&=\, \imath\,(a_m-a_{-m})\\ \end{split}. $$

The concept Fourier series can be expanded to a base interval of a period T instead of $2\pi$ by substituting $x$ with $x = \frac{2\pi}{T}t$.

$$ g_{\rm F}(t) = f_{\rm F}(\frac{2\pi}{T}t) \,=\, \frac{1}{2}c_0+\sum_{m = 1}^{\infty}c_m \,\cos(m\frac{2\pi}{T}t)+\sum_{m = 1}^{\infty}s_m \,\sin(m\frac{2\pi}{T}t) $$

where $$ c_0 \,=\,\frac{1}{\pi}\int_{-\pi}^{\pi}f(\frac{2\pi}{T}t)\,dx \,=\, \frac{2}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}}g(t)\,dt. \\ c_m \,=\,\frac{1}{\pi}\int_{-\pi}^{\pi}f(\frac{2\pi}{T}t)\,\cos(m\frac{2\pi}{T}t)\,dx \,=\, \frac{2}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}}g(t)\,\cos(m\frac{2\pi}{T}t)\,dt \qquad m \in \mathbb{N_0}\\ s_m \,=\,\frac{1}{\pi}\int_{-\pi}^{\pi}f(\frac{2\pi}{T}t)\,\sin(m\frac{2\pi}{T}t)\,dx \,=\, \frac{2}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}}g(t)\,\sin(m\frac{2\pi}{T}t)\,dt \qquad m \in \mathbb{N_0}\\ $$

or

$$ g_{\rm IF}(t) = f_{\rm IF}(\frac{2\pi}{T}t) \,=\, \sum_{m = -\infty}^{\infty}a_m \,e^{\imath m\frac{2\pi}{T}t} $$

$$ a_m \,=\, \frac{1}{2\pi}\int_{-\pi}^{\pi}f(\frac{2\pi}{T}t)e^{-\imath m\frac{2\pi}{T}t}\,dx\,=\,\frac{1}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}}g(t)e^{-\imath m\frac{2\pi}{T}t}\,dt \qquad\forall m\in\mathbb{Z}.$$

The series again converges under the same criteria as before and the relations between the coefficients of the complex or real Fourier coefficients from equation equation stay valid.

One nice example is the complex, scaled Fourier series of the scaled shah function ➞ $III_{T^{-1}}(x)\,=III\left(\frac{x}{T}\right)\,=\sum_{l=-\infty}^{+\infty} T \delta\left(x-l T\right)$. Obviously, the period of this function is $T$. The Fourier coefficients (matched to a period of $T$) is calculated as

$$ \begin{split} a_m \,&= \,\frac{1}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}}\left(\sum_{l=-\infty}^{+\infty}T\delta\left(x-l T\right)\right)e^{-\imath m \frac{2\pi}{T} x}\,dx\\ &=\,\frac{1}{T} \int_{-\frac{T}{2}}^{\frac{T}{2}} T \delta\left(x\right)e^{-\imath m \frac{2\pi}{T} x}\,dx\\ &=\,1 \end{split} \forall m\in\mathbb{Z}.$$

It follows that

$$ \begin{split} III_{T^{-1}}(x)\,=III\left(\frac{x}{T}\right)\,=\,\sum_{m=-\infty}^{+\infty} e^{\imath m\frac{2\pi}{T} x t} \end{split} .$$

2.3.2 Example

Next we will quickly demonstrate how to decompose a signal into its Fourier series. An easy way to implement this numerically is to use the trapezoidal rule to approximate the integral. Thus we start by defining a function that computes the coefficients of the Fourier series using the complex definition


In [ ]:
def FS_coeffs(x, m, func, T=2.0*np.pi):
    """
    Computes Fourier series (FS) coeffs of func
    Input: 
        x = input vector at which to evaluate func
        m = the order of the coefficient
        func = the function to find the FS of
        T = the period of func (defaults to 2 pi)
    """
    # Evaluate the integrand
    am_int = func(x)*np.exp(-1j*2.0*m*np.pi*x/T)
    # Use trapezoidal integration to get the coefficient
    am = np.trapz(am_int,x)
    return am/T

That should be good enough for our purposes here. Next we create a function to sum the Fourier series.


In [ ]:
def FS_sum(x, m, func, period=None):
    # If no period is specified use entire domain
    if period is None:
        period = np.abs(x.max() - x.min())
    
    # Evaluate the coefficients and sum the series
    f_F = np.zeros(x.size, dtype=np.complex128)
    for i in xrange(-m,m+1):
        am = FS_coeffs(x, i, func, T=period)
        f_F += am*np.exp(2.0j*np.pi*i*x/period) 
    return f_F

Let's see what happens if we decompose a square wave.


In [ ]:
# define square wave function
def square_wave(x):
    I = np.argwhere(np.abs(x) <= 0.5)
    tmp = np.zeros(x.size)
    tmp[I] = 1.0
    return tmp

# Set domain and compute square wave
N = 250
x = np.linspace(-1.0,1.0,N)

# Compute the FS up to order m
m = 10
sw_F = FS_sum(x, m, square_wave, period=2.0)

# Plot result
plt.figure(figsize=(15,5))
plt.plot(x, sw_F.real, 'g', label=r'$ Fourier \ series $')
plt.plot(x, square_wave(x), 'b', label=r'$ Square \ wave $')
plt.title(r"$FS \ decomp \ of \ square \ wave$",fontsize=20)
plt.xlabel(r'$x$',fontsize=18)
plt.ylim(-0.05,1.5)
plt.legend()

Figure 2.8.1:

As can be seen from the figure, the Fourier series approximates the square wave. However at such a low order (i.e. $m = 10$) it doesn't do a very good job. Actually an infinite number of Fourier series coefficients are required to fully capture a square wave. Below is an interactive demonstration that allows you to vary the parameters on the Fourier series decomposition. Note in particular what happens if we make the period too small. Also feel free to apply it to functions other than the square wave (but make sure to adjust the domain accordingly.


In [ ]:
def inter_FS(x,m,func,T):
    f_F = FS_sum(x, m, func, period=T)
    plt.plot(x,f_F.real,'b')
    plt.plot(x,func(x),'g')
    

interact(lambda m,T:inter_FS(x=np.linspace(-1.0,1.0,N),m=m,func=square_wave,T=T),
                m=(5,100,1),T=(0,2*np.pi,0.5)) and None

Figure 2.8.2: