# Fourier Series representation using Python
Fourier series, in my view, is among the most revolutionary ideas in mathematics. It was [Joseph Fourier](https://en.wikipedia.org/wiki/Joseph_Fourier)'s imaginative brilliance, that he managed to show the world how a periodic discontinuous signal like a square wave can be represented using something as continous as sinusoidal waves. And that too all using his mathematics skills and imagination, when modern visualization tools weren't available.
In this post I'm going to visualize the Fourier series representation of a simple square wave using two of the most common Python libraries: **NumPy** and **Matplotlib**.
But before we jump onto the visualization part, we need to know what does the Fourier series representation of a periodic square wave look like mathematically. This is important so as to know exactly how (and why), the Fourier Series will behave.
## Mathematical Representation
Let us consider the square wave:
The time period of this simple square wave is T=4.
The general expression for a trigonometric Fourier series, looks something like this$^{[1]}$ :
$$x(t) = a_{0} + \sum_{k=1}^{\infty} [a_{k}cos(k\omega_{o}t)+b_{k}sin(k\omega_{o}t)] ...(1)$$
The Fourier coefficients in the above expression are given as:
$$a_{0} = \frac{1}{T}\int_{0}^{T}x(t)dt$$
$$a_{k} = \frac{2}{T}\int_{0}^{T}x(t)cos(k\omega_{o}t)dt$$
$$b_{k} = \frac{2}{T}\int_{0}^{T}x(t)sin(k\omega_{o}t)dt$$
By applying these formulas for the square wave above, we get the fourier coefficients as:
$$a_{0} = \frac{1}{2}$$
$$a_{k} = \frac{2sin(\frac{k\pi}{2})}{k\pi}$$
$$b_{k} = 0$$
I've puposedly not derived these expressions here(which can be easily found in standard textbooks$^{[2]}$) so as to keep the post more oriented towards the visualization part.
Note that the value $b_{k}=0$ indicates that there will be no sine terms in the Fourier representation of this signal which is essentially because this signal is an [even signal](https://en.wikipedia.org/wiki/Even_and_odd_functions).
Putting different values of $k$, we get the Fourier coefficients of the cosine terms as:
$a_{1}=\frac{2}{\pi},a_{2}=0, a_{3}=\frac{-2}{3\pi},a_{4}=0,a_{5}=\frac{2}{5\pi},a_{6}=0$ ...and so on.
Note, that the even terms (except $a_{0}$) are all $0$.
So finally, plugging these values in the equation (1) above, we get the tigonometric Fourier series representation of the square wave as:
$$x(t)=\frac{1}{2}+\frac{2cos(\omega_{o}t)}{\pi}+\frac{-2cos(3\omega_{o}t)}{3\pi}+\frac{2cos(5\omega_{o}t)}{5\pi}+...$$
Note, that there are no sine terms as all the $b_{k}=0$.
So, yeah, after all this heavy mathematics, we've achieved the final Fourier represenation of our square wave. But before, we can jump to visualization we need to see a few more points:
* First of all see that the term $a_{0}=\frac{1}{2}$ actually represents the DC component of the signal $x(t)$.
Next, by observation one can see that a general term of the Fourier series calculated above will look something like this:
$$X_{k}={(-1)}^{\frac{k-1}{2}}\frac{2cos(\frac{2\pi kt}{T})}{k\pi}...(2)$$
* As we put different values of $k$ we get different cosine terms.
Like when $k=1$, $X_{1}=\frac{2cos(\frac{2\pi t}{T})}{\pi}$
when $k=3$, $X_{3}=-\frac{2cos(\frac{2\pi 3t}{T})}{3\pi}$
when $k=5$, $X_{5}=\frac{2cos(\frac{2\pi 5t}{T})}{5\pi}$
and so on...
* These cosine terms are called *harmonically realted sinusoids*. The term 'harmonically-realted' arises from the fact that the frequency of each of these sinusoids is actually a multiple of the fundamental frequency $\omega_{o}$.
* So, this concludes the **basic idea of Fourier series representation**: A periodic signal can be represented using sum of harmonically related sinusoids$^{[3]}$.
For the sake of convinience, in this post, 'harmonically-related sinusoids' will be simply reffered to as 'harmonics'.
We now turn to our main task: we need to see through *actual visualization* whether the mathematical representation obtained above for our sqaure wave, *really represents* the square wave? More precisely, we *wish to see how much* does this representation resembles the actual square wave.
## Visualization
Theoretically, the idea is that, as more and more terms (harmonics) get added to the Fourier representation, the more closely the function should represent the original square wave.
To visualize this idea, we'll plot the mathematical function obtained above using Python and add more and more harmonics to see how closely it resembles our actual square wave.
Primarily we'll be using NumPy and Matplotlib- two of the most popular Python libraries. You can install these libraries to run it locally on your system. You can install [Anaconda](https://www.anaconda.com/), which is a python distribuition and comes with built-in python packages and a very powerful tool- Jupyter notebooks. If you wish to save yourselves from all of this installation, simply head on to colab.research.google.com and start writing your code!
---
Look at the code below for reference. The first three lines simply import the required libraries. The 4th line initializes the time axis `t`. We can adjust the range and scale of $t$ according to our will. Next we set the time period `T=4`. Next we initialize the DC component (calculated as $a_{0}=\frac{1}{2}$). Next we set the `amplitude` variable which basically represnts the value of the signal at different time instants. The $k$ in the 8th line stands for the same $k$ as we've referred to before- the harmonic number. Now comes the main part. How do we add more and more harmonics to our `amplitude`? We do this by using a while loop which starts from line 10. The value of `i` will basically indicate the number of harmonics we're adding. So, if we want to add n number of harmonics we'll iterate n times through this loop. Inside this loop at line number 11, we have initialized the harmonic term `X` (the expression which has already been mentioned in equation (2)). As we iterate over this loop the value of `amplitude` gets updated. Finally in Line 15 we plot the graph of this function using the `plot` function of Matplotlib.
```python=
import numpy as np
import matplotlib.pyplot as plt
from math import pi as p
t=np.arange(-5,5,0.0001) #time axis
T=4 # Time Period
X0=0.5 # DC Component
amplitude=X0
k=1
i=0
while i<3: #Change the value of i to change the number of harmonics added
X=(np.power(-1, (k-1)/2))*(2/(p*k))*np.cos((k*2*p*t)/T)
amplitude=amplitude+X
k+=2
i+=1
plt.plot(t,amplitude,color='red')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.title('Fourier Series representation of Periodic Signal')
plt.grid()
plt.show
```
Executing the above code snippet will give you the following graph:
This wave above represents the Fourier series representation of square wave (Eqn 2).
Okay, so if you're honest to yourselves, you'll instantly realise that this waveform above actually resembles a bit to the square wave. In this example we've added just three harmonics (`i<3`) and the Fourier series waveform has almost taken the shape of the sqaure wave.
However, it still doesn't completely match the square wave; so let's increase the number of harmonics by increasing the number of times the loop is run: `i<7`
That's definitely closer to the square wave as compared to the previous figure. Let's increase the harmonics further:`i<21`
That's even closer. So we can observe that as we're increasing the number of harmonics added to the Fourier series, the wavefom is becoming sharper and closer to the actual square wave.
You can tweak the value of parameter in `i<n` as many times as you wish to change the number of harmonics and see for yourselves how the shape of the Fourier representation changes.
### Observations
By observing these visualizations we can make some inferences:
* As we're increasing the number of times the while loop will run, we're adding more and more harmonics and more closely the graph is reprsenting our actual square wave.
* It can be observed from the first couple of figures that the lower harmonics have almost taken the shape of the square wave and adding higher harmonics just increases the sharpness of this shape.
* Notice that at the point of discontinuities in the waveform, we can observe some overshoots. This is actually a property of Fourier series known as the [Gibbs Phenomenon](https://en.wikipedia.org/wiki/Gibbs_phenomenon).
## Footnotes
[1] For Derivation see: https://lpsa.swarthmore.edu/Fourier/Series/DerFS.html
[2] See Example 3.5, Page no. 195 of _Signals and Systems_ by Alan V. Oppenheim.
[3] In more mathematical form the word 'sum' is replaced by the term 'linear combination'
Sign in with Wallet
Connect another wallet