In this article, we will walk through the origins of the Fourier transform: the **Fourier Series**. The Fourier series takes a periodic signal $x(t)$ and describes it as a sum of sine and cosine waves. Noting that sine and cosine are themselves periodic functions, it becomes clear that $x(t)$ is also a periodic function.

Mathematically, the Fourier series is described as follows. Let $x(t)$ be a periodic function with period $T$, i.e.

$$x(t)=x(t+nT), n\in\mathbb{Z}.$$Then, we can write $x(t)$ as a Fourier series by

$$x(t)=\frac{a_0}{2}+\sum_{n=1}^{\infty}a_n\cos(2\pi \frac{nt}{T})+b_n\sin(2\pi\frac{nt}{T}),$$where $a_n$ and $b_n$ are the coefficients of the Fourier series. They can be calculated by $$\begin{align}a_n&=\frac{2}{T}\int_0^Tx(t)\cos(2\pi \frac{nt}{T})dt\\ b_n&=\frac{2}{T}\int_0^Tx(t)\sin(2\pi \frac{nt}{T})dt\end{align}.$$

Note that for a function with period $T$, the frequencies of the sines and cosines are $\frac{1}{T}, \frac{2}{T}, \frac{3}{T}, \dots$, i.e. they are multiples of the fundamental frequency $\frac{1}{T}$, which is the inverse period duration of the function. Therefore the frequency $\frac{n}{T}$ is called the $n$th *harmonic*. The name *harmonic* stems from the fact for the human ear frequencies with integer ratios sound "nice", and the frequencies are all integer multiples of the fundamental frequency.

Let us verify the calculation of the Fourier coefficients and the function reconstruction numerically. First, we define some functions with period $T=1$ that we want to expand into a Fourier series:

```
Fs = 10000
func1 = lambda t: (abs((t%1)-0.25) < 0.25).astype(float) - (abs((t%1)-0.75) < 0.25).astype(float)
func2 = lambda t: t % 1
func3 = lambda t: (abs((t%1)-0.5) < 0.25).astype(float) + 8*(abs((t%1)-0.5)) * (abs((t%1)-0.5)<0.25)
func4 = lambda t: ((t%1)-0.5)**2
t = np.arange(-1.5, 2, 1/Fs)
plt.subplot(221); plt.plot(t, func1(t))
plt.subplot(222); plt.plot(t, func2(t))
plt.subplot(223); plt.plot(t, func3(t))
plt.subplot(224); plt.plot(t, func4(t))
```

`fourierSeries`

that performs the calculation of the Fourier series coefficients:

```
def fourierSeries(period, N):
"""Calculate the Fourier series coefficients up to the Nth harmonic"""
result = []
T = len(period)
t = np.arange(T)
for n in range(N+1):
an = 2/T*(period * np.cos(2*np.pi*n*t/T)).sum()
bn = 2/T*(period * np.sin(2*np.pi*n*t/T)).sum()
result.append((an, bn))
return np.array(result)
```

And use it to calculate the coefficients up to the 20th order for the first function:

```
t_period = np.arange(0, 1, 1/Fs)
F = fourierSeries(func1(t_period), 20)
plt.subplot(121); plt.stem(F[:,0])
plt.subplot(122); plt.stem(F[:,1])
```

We make two observations here: First, we see that $a_n=0$. From the plot for function 1 we see that it is an odd function, i.e. $x(t)=-x(-t)$. In this case, the Fourier series only contains odd functions, which are solely the terms including the sines (since $\sin(x)=-\sin(-x)$). Second, the Fourier coefficients $b_n$ decay slowly with speed $1/n$ and every 2nd $b_n$ is zero. We can explain this by relating it to knowledge of the Fourier Transform: The function is a sum of two rectangular functions of width $\frac{1}{2}s$. We know the Fourier transform of such a rectangle is a sinc-function, which has zeros at a distance of $2Hz$. Furthermore, the magnitude of the sinc-function decays with $1/f$. This is very in line with the obtained coefficients: They decay with $1/n$ and every 2nd values is zero.

Let us now look at the reconstruction of the signal, i.e. we calculate $x(t)$ from its Fourier series coefficients up to a given order. We write a function `reconstruct`

to do this for us:

```
def reconstruct(P, anbn):
result = 0
t = np.arange(P)
for n, (a, b) in enumerate(anbn):
if n == 0:
a = a/2
result = result + a*np.cos(2*np.pi*n*t/P) + b * np.sin(2*np.pi*n*t/P)
return result
```

Let's have a look at the reconstructed signal for our rectangular function up to the 20th harmonic:

```
F = fourierSeries(func1(t_period), 100)
plt.plot(t_period, func1(t_period), label='Original', lw=2)
plt.plot(t_period, reconstruct(len(t_period), F[:20,:]), label='Reconstructed with 20 Harmonics');
plt.plot(t_period, reconstruct(len(t_period), F[:100,:]), label='Reconstructed with 100 Harmonics');
```

As we can see, the reconstructed signal roughly follows the original, and the more harmonics are used, the better. However, we also see that especially in the region of the jump at $t=0.5$, the reconstructed signal is not exact. Instead, the reconstructed signal significantly fluctuates at this position. This phenomenom is called Gibbs Phenonenom and describes the fact that the Fourier series has large oscillations around jump discontinuities. In particular, the height of overshooting or undershooting does not depend on the number of harmonics and is roughly 9% of the jump height. However, the duration of the oscillations decreases with the number of harmonics, eventually leading to a correct approximation in the limit for infinitely many harmonics.

Let us now have a look at the Fourier Series of some functions, and how their approximation by the Fourier series appears for different number of Harmonics:

```
def showHarmonics(period, N):
"""Calculate the Fourier Series up to N harmonics, and show the reconstructed signal."""
F = fourierSeries(period, N+1)
plt.subplot(231); plt.stem(F[:,0])
plt.subplot(234); plt.stem(F[:,1])
plt.subplot(132)
T = len(period)
t = np.arange(T)/T
result = 0
for n, (an, bn) in enumerate(F):
if n == 0:
an = an/2
cos_part = an*np.cos(2*np.pi*n*t)
sin_part = bn*np.sin(2*np.pi*n*t)
plt.plot(t, cos_part)
plt.plot(t, sin_part)
result = result + cos_part + sin_part
plt.subplot(133)
t2 = np.arange(2*T)/T
plt.plot(t2, np.tile(period, 2))
plt.plot(t2, np.tile(result, 2))
```

```
```

```
```

```
```

```
```

- The Fourier Series decomposes a periodic function with period T into sines and cosines with frequencies $\frac{n}{T}, n=0,1,2,\ldots$ which are called the $n$th harmonics of the signal. The more harmonics are used, the more accurate can a function be described.
- For
even functions, i.e. $x(t)=x(-t)$, the Fourier seriesonly consists of cosines. Forodd functions, i.e. $x(t)=-x(-t)$, the Fourier seriesonly consists of sines.- At jump discontinuities, the value of the Fourier series is in the middle of the jump. Around the jump, overshooting of the series occurs, which is called
Gibbs Phenonemon. The amount of overshooting does not reduce with the number of harmonics, but the duration reduces.- Smooth functions need less harmonics to be accurately described by the Fourier series.

In a future article, we will investigate the relation of the Fourier series to the Fourier transform, and analyze sounds of instruments according to their structure of the harmonics. Subscribe to our newsletter on the right to not miss upcoming posts!

Do you have questions or comments? Let's dicuss below!