January 11, 2024
Any periodic function can be expressed as an infinite series of sine and cosine functions. Specifically, a function $f(x)$ that is periodic on the interval $-\pi<x<\pi$ can be expressed as: $$ \begin{aligned} f(x) &= \frac{a_0}{2} + a_1\cos x + a_2\cos 2x + a_3\cos 3x + \dots + b_1\sin x + b_2\sin 2x + b_3\sin 3x + \dots\\ &=\frac{a_0}{2} + \sum_{n=1}^\infty\left(a_n\cos nx + b_n\sin nx\right) \end{aligned} $$ The coefficients are determined by the following integrals: $$ \begin{aligned} a_0 &= \frac{1}{\pi}\int_{-\pi}^\pi f(x)\, dx\\ a_n &= \frac{1}{\pi}\int_{-\pi}^\pi f(x)\cos nx\, dx\\ b_n &= \frac{1}{\pi}\int_{-\pi}^\pi f(x)\sin nx\, dx \end{aligned} $$ where $n = 1, 2, 3, \dots$
Below, we attmept to show how a superposition of sinusoidal functions can be used to construct a square wave. Recall that a square wave of ampltude 1 and period $2\pi$ can be expressed as: $$ \begin{aligned} Sq(x) &= \frac{1}{2} + \frac{2}{\pi}\left(\sin x + \frac{1}{3}\sin 3x + \frac{1}{5}\sin 5x + \frac{1}{7}\sin 7x + \dots\right)\\ &= \frac{1}{2} + \sum_{n,~\mathrm{odd}}\frac{\sin nx}{n} \end{aligned} $$ The function $Sq(x)$ above is for a square wave that oscillates between 0 and 1.
Import a couple of required modules.
import numpy as np
import matplotlib.pyplot as plt
Define a set of $x$ values that range from $-2\pi$ to $+2\pi$. This range of values will show two full periods of the square wave.
xx = np.linspace(-2*np.pi, 2*np.pi, 1000)
First, we show only the first two terms in the Fourier series (the $a_0$ and $b_1$ terms). Of course, this results in a pure sine wave with a constant offset.
fcn = 1/2 + (2/np.pi)*np.sin(xx)
fig, ax = plt.subplots()
ax.plot(xx, fcn, 'r')
plt.xlabel('x')
plt.ylabel('fcn');
Next, we add the $b_3$ term. The first plot is the square wave approximation using the first three non-zero terms of the Fourier expansion.
fcn = 1/2 + (2/np.pi)*(np.sin(xx) + (1/3)*np.sin(3*xx))
plt.plot(xx, fcn, 'b')
plt.xlabel('x')
plt.ylabel('fcn');
This plot overlays the red and blue plot above in order to highlight the evolution towards a better and better approximation of a square wave.
ax.plot(xx, fcn, 'b')
fig
We now include the $a_0$, $b_1$, $b_3$, and $b_5$ terms.
fcn = 1/2 + (2/np.pi)*(np.sin(xx) + (1/3)*np.sin(3*xx) + (1/5)*np.sin(5*xx))
plt.plot(xx, fcn, 'g')
plt.xlabel('x')
plt.ylabel('fcn');
ax.plot(xx, fcn, 'g')
fig
We now include the $a_0$, $b_1$, $b_3$, $b_5$, and $b_7$ terms.
fcn = 1/2 + (2/np.pi)*(np.sin(xx) + (1/3)*np.sin(3*xx) + (1/5)*np.sin(5*xx) + (1/7)*np.sin(7*xx))
plt.plot(xx, fcn, 'orange')
plt.xlabel('x')
plt.ylabel('fcn');
ax.plot(xx, fcn, 'orange')
fig
We'll now use a for loop to include all non-zero terms in the Fourier expansion up to $n = 29$.
fcn = 1/2
for n in range(1, 30, 2):
fcn = fcn + (2/np.pi)*(1/n)*np.sin(n*xx)
As the plot below shows, we are now getting closer to a square wave shape.
plt.plot(xx, fcn, 'cyan')
plt.xlabel('x')
plt.ylabel('fcn');
ax.plot(xx, fcn, 'cyan')
fig
Finally, we use the same for loop, but now use it to construct a Fourier series that includes a non-zero terms up to $n = 9999$. The result is a nearly perfect square wave.
fcn = 1/2
for n in range(1, int(1e4), 2):
fcn = fcn + (2/np.pi)*(1/n)*np.sin(n*xx)
plt.plot(xx, fcn, 'k')
plt.xlabel('x')
plt.ylabel('fcn');