The Fourier Series

Any periodic function can be expressed as a sum of sines and cosines. Period.

That's basically what the Fourier Series is. All thanks to Joseph Fourier's study of Heat and it's propagation over matter. Fourier realized that temperature in materials could be expressed as a combination of sine and cosine waves.

"Heat, like gravity, penetrates every substance of the universe, its rays occupy all parts of space. The object of our work is to set forth the mathematical laws which this element obeys." ~ Joseph Fourier.

Definition

The Fourier Series allows any periodic function f(x) to be expressed as a (possibly infinite) sum of sines and cosines, as follows:

Eq 1. \[ f(x) = \frac{a_0}{2} + \sum_{n=1}^\infty \, [a_n \cos(nx) + b_n \sin(nx)] \]

For this periodic function f(x), integrable in the interval [-π, π], the numbers a_n and b_n are known as the Fourier Coefficients of f, and defined as:

Eq 2. \[ a_n = \frac{1}{\pi}\int_{-\pi}^\pi f(x) \cos(nx)\, dx, \quad n \ge 0 \]

Eq 3. \[ b_n = \frac{1}{\pi}\int_{-\pi}^\pi f(x) \sin(nx)\, dx, \quad n \ge 1 \]

Application of Fourier Series in Python

The Fourier series allows a periodic function to be expressed in terms of sines and cosines. This sum is often infinite and Fourier coefficients need to be calculated too. This process is repetitive and normally a manual calculation is not practical.

Python allows to automate these calculations and to obtain a graphic result of each approximation. These are called approximations because, in many cases, the periodic function will be equal to an infinite sum of sines and cosines. Even though an infinite calculation is not always possible, an approximation can be made with a good number of iterations (100, 500, etc).

Firstly, we'll need the NumPy, SciPy, and matplotlib libraries. It's worth highlighting that the integrate module from SciPy allows to perform numerical integrations.

import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt
import scipy as sp
import matplotlib as mpl

Now we proceed to calculate the first sum of equation 2. To do this, the value of a0 is calculated from equation 1.

 a0, err0 = integrate.quad(lambda x: x,-sp.pi,sp.pi)

Here, the quad function allows us to integrate a function between two points. In this example, the function f(x) is equal to x, that is, a straight line with period 2π:

\[ f(x) = x, \quad \mathrm{for } -\pi < x < \pi, \]

\[ f(x + 2\pi) = f(x), \quad \mathrm{for } -\infty < x < \infty. \]

The lower limit is -π and the upper one is π according to equation 2. The 'quad' function returns the result of the integral and an estimate of the absolute error of the integration; These results are stored in the variables a0 and err0 respectively. We'll then divide the integral by π, according to equation 2.

a0 = (1/sp.pi)*a0
y = a0*0.5
z = np.arange(0,5*np.pi,0.01)

To graph the approximations to the Fourier series it is necessary to tabulate the values of x and f(x). These approximations have the form

Eq. 4 \[ f(x)\approx\frac{a_0}{2} + \sum_{n=1}^N \, [a_n \cos(nx) + b_n \sin(nx)], \quad N \ge 0. \]

Unlike equation 1, equation 4 is simply an approximation to the Fourier series. By varying the value of N, increasingly accurate approximations are obtained. When N tends to infinity, the Fourier series of the function f(x) is obtained.

"y" is a vector that will store all the values (the ordinates) of f(x). The variable y is initialized with the value a_0/2 because, as can be seen in equation 4 (and 3) a_0/2 is added to all the values of the ordinates.

On the other hand, we'll the z vector to store the values of the abscissas. It will store values from 0 to 5 semi-periods (5p), with steps of 0.01, as follows:

z=[0,0.01,0.02,...,5π]

The range [0,5p] will allow us to graph only 5 semi-periods of the series.

We'll the proceed to calculate equation 4 will be added to the y array. To do this, the values of the Fourier constants (equations 2 and 3) are calculated.

As seen in the equations to find the Fourier coefficients a_n and b_n, the integration function must be the function f(x) multiplied by cos(nx) and sin(nx) respectively. For the case of a_0 and b_0 the cosine and sine are 1 and 0, respectively.

"y += ai*np.cos(i*z) + bi*np.sin(i*z)" updates all the entries of vector y (which, initially, are all equal to a_0/2), with a new pair or linear combination of sines and cosines. It's all about adding adding sines and cosines.

\[ z=\left [\begin{matrix} 0 \\ 0.01 \\ 0.02 \\ \vdots \\ 5\pi \\ \end{matrix}\right ] \ \wedge \ y=\left [\begin{matrix} \frac{a_0}{2} + a_1 \cos(1\cdot0) + b_1 \sin(1\cdot0) + \cdots + a_i \cos(i\cdot0) + b_i \sin(i\cdot0) + \cdots + a_N \cos(N\cdot0) + b_N \sin(N\cdot0) \\ \frac{a_0}{2} + a_1 \cos(1\cdot0.01) + b_1 \sin(1\cdot0.01) + \cdots + a_i \cos(i\cdot0.01) + b_i \sin(i\cdot0.01) + \cdots + a_N \cos(N\cdot0.01) + b_N \sin(N\cdot0.01) \\ \frac{a_0}{2} + a_1 \cos(1\cdot0.02) + b_1 \sin(1\cdot0.02) + \cdots + a_i \cos(i\cdot0.02) + b_i \sin(i\cdot0.02) + \cdots + a_N \cos(N\cdot0.02) + b_N \sin(N\cdot0.02) \\ \vdots \\ \frac{a_0}{2} + a_1 \cos(1\cdot5\pi) + b_1 \sin(1\cdot5\pi) + \cdots + a_i \cos(i\cdot5\pi) + b_i \sin(i\cdot5\pi) + \cdots + a_N \cos(N\cdot5\pi) + b_N \sin(N\cdot5\pi) \\ \end{matrix}\right ] \]

Finally, the abscissas (vector z) and ordinates (vector y) are plotted.

import numpy as np
import scipy as sp
from scipy import integrate
import matplotlib.pyplot as plt
import matplotlib.font_manager as fm

# Define fonts to use on Matplotlib plots.
TitleFont = fm.FontProperties(fname='Quicksand-Regular.otf', size=25)

# Calculation of coefficients
a0, err0 = integrate.quad(lambda x: np.exp(x), -sp.pi, sp.pi)
a0 = (1/sp.pi)*a0
y = a0*0.5
z = np.arange(0, 5*np.pi, 0.01)

# Create figure
fig = plt.figure(figsize=(5, 7))

# Number of subplots
N = 140
iterations = [1, 2, 5, 25, 120]

# Iterate over subplots
for i in range(1, N + 1):
    ai, erri = integrate.quad(lambda x: x * sp.cos(i*x), -sp.pi, sp.pi)
    ai = (1/sp.pi)*ai
    bi, erri = integrate.quad(lambda x: x * sp.sin(i*x), -sp.pi, sp.pi)
    bi = (1/sp.pi)*bi
    y += ai*np.cos(i*z) + bi*np.sin(i*z)
    
    if i in iterations:
      num_of_graphs = len(iterations)
      index = iterations.index(i) + 1

      # Add subplot
      ax = fig.add_subplot(num_of_graphs, 1, index)
      
      # Remove x-axis and y-axis numbers
      ax.xaxis.set_visible(False)
      ax.yaxis.set_visible(False)
      
      # Remove the frame around the subplot
      ax.set_frame_on(False)

      # Set title for each subplot
      ax.set_title(f'Iteration #{i}', fontproperties=TitleFont, fontsize=15)

      # Plot data
      color = '#f04859'
      ax.plot(z, y, color)

plt.subplots_adjust(wspace=2, hspace=0.5)
plt.show()

Each time the code is executed, the N-th approximation will be graphed. To choose which approximation you want to graph (or how many iterations you want to do), you just have to choose the N value in the line "N = a".