# 13. Fourier series: the “bones” of a function

There is often more to a plot than immediately meets the eye, and in the life of a scientist one of the great joys is to glean more information than what is on the surface.

Here we will show one of my favorite tools for extracting further
information from data: the *Fourier series*. The result of this
transformation is referred to as the *Fourier spectrum* of a signal.

Let us immediately dive in with a surprising example: breaking the
*square wave* down into a sum of very round sin waves.

## 13.1. Fourier analysis: the square wave

*Square waves* are used as “clocks” in electronics: the edges can
trigger logic circuits to do an action simultaneously with other
circuits.

Let us start gnuplot and plot a \(\sin\) wave, then look at some \(\sin\) waves with higher frequencies. Then add them together and see what you get:

```
$ gnuplot
# and at the gnuplot> prompt:
set samples 1000
plot [] [-1.2:1.2] sgn(sin(x))
replot (4.0/pi) * sin(x)
replot (4.0/pi) * (1.0/3)*sin(3*x)
replot (4.0/pi) * (1.0/5)*sin(5*x)
replot (4.0/pi) * (sin(x) + (1.0/3)*sin(3*x) + (1.0/5)*sin(5*x))
## now look at the summed-up plot by itself:
plot (4.0/pi) * (sin(x) + (1.0/3)*sin(3*x) + (1.0/5)*sin(5*x)) t '5x'
```

You should start seeing that you go from a \(\sin\) wave to something that looks a bit more like a square wave. shows what the individual \(\sin\) waves look like and shows how you can add up to 13 of them and get something that starts to look quite square instead of wavy.

This kind of series that sums up to the square function is called the
*Fourier series* for the square wave. The full expression is:

where \(\omega = 2\pi f\) is the *angular frequency*.

Note

Here the instructor should spend some time talking about the angular frequency. The nature of this conversation depends on how much exposure to physics the students have had. The simplest discussion of it could be that instead of counting how many times a periodic function repeats in a unit of time, we count how much angle a circular movement covers in a unit of time. Of course the angle is measured in radians, which is how we get \(\omega = 2\pi f\).

More generally *you can approximate any periodic function* with period
\(T\) with the series:

where \(\omega = \frac{2\pi}{T}\).

Terminology:

- fourier series
The way of writing the function broken down as a sum of sin and cos waves.

- fourier coefficients
The values \(a_k \; {\rm and} \; b_k\) are the

*fourier coefficients*for this function.- signal
One of the most active areas of science and engineering in which Fourier analysis is used is called

*signal processing*, where we record electrical signals and analyze them in various ways. In this kind of setting we often use the word*signal*to mean the function we are discussing.- fourier spectrum
The ensemble of all of the fourier coefficients, \({a_k, b_k}\), is called the

*fourier spectrum*of the signal.

## 13.2. Derivation of the Fourier series for periodic signals

We will talk here about periodic functions of a single real variable (we’ll call it \(t\)) with period \(T\). In these functions \(f(t + T) = f(t)\).

For this arbitrary periodic function \(f(t)\), the Fourier series is given by the expression shown above:

and the coefficients are given by:

How do you derive these expressions for the \(a_n \; {\rm and} \; b_n\)? Start by remembering this fact from integration:

(Similar formulas apply to the \(\sin()\) function.)

This allows us to see that multiplying our series by
\(\cos(m\omega t)\) will cause some simplifications when \(m
\ne n\). We won’t look at the full derivation, which can be found as a
straightforward result of searching for *fourier series derivation* -
I just wanted to point out that this kind of integral over a full
period of trigonometric functions allows us to pick out the
coefficients.

Note

The fourier series is also sometimes written as:

where the coefficients are \(A_n \; {\rm and} \; \phi_n\). The trigonometric identity for \(\cos(\alpha + \beta)\) will quickly show you that you can translate this into our sin + cos formula, mapping these coefficients to the \(a_n \; {\rm and} \; b_n\).

Caution

The Fourier transform is so often used in science, mathematics, engineering, and music that scholars have come up with many different conventions for how to write it. All these different conventions have slightly different notations or choices on whether to put certain factors in certain places, but they all end up having the same information. This means that when you read a paper where Fourier analysis is used, you have to look carefully to understand what convention they use. If the author did not specify a convention, then you should resent their sloppiness.

## 13.3. A tour of functions and their Fourier series

We now apply our formulas for \(a_n \; {\rm and} \; b_n\) to a few common types of periodic signals. Remember that if it is periodic then we can just act like the function is defined on the interval \([0:2\pi]\) or \([-\pi:\pi]\). To get the full extent of the function we just duplicate that single period forever in both directions.

### 13.3.1. Sawtooth wave

The coefficients are:

So the series is:

Let us plot this (using a period of \(T = 1\)) with:

```
reset
set grid
set samples 1000
set xrange [-2:2]
set yrange [-1.3:1.3]
# set terminal qt lw 2
plot 0 lw 1 lc 'black'
replot [-1:1] 1./2 + x/2 lw 2 lc 'green'
set arrow from -1, 0 to -1, 1 nohead lc 'green' lw 2
replot [1:2] 1./2 + x/2 - 1 lw 2 lc 'green'
set arrow from 1, 0 to 1, 1 nohead lc 'green' lw 2
replot [-2:-1] 1./2 + x/2 + 1 lw 2 lc 'green'
# that was the sawtooth wave. pause to contemplate it, and then
# start with the Fourier series:
replot 1./2
replot 1./2 - 1/(pi)*( (-1)**1 * sin(1*pi*x)/1)
replot 1./2 - 1/(pi)*( (-1)**1 * sin(1*pi*x)/1 + (-1)**2 * sin(2*pi*x)/2)
replot 1./2 - 1/(pi)*( (-1)**1 * sin(1*pi*x)/1 + (-1)**2 * sin(2*pi*x)/2 + (-1)**3 * sin(3*pi*x)/3 + (-1)**4 * sin(4*pi*x)/4 + (-1)**5 * sin(5*pi*x)/5 + (-1)**6 * sin(6*pi*x)/6)
replot 1./2 - 1/(pi)*( (-1)**1 * sin(1*pi*x)/1 + (-1)**2 * sin(2*pi*x)/2 + (-1)**3 * sin(3*pi*x)/3 + (-1)**4 * sin(4*pi*x)/4 + (-1)**5 * sin(5*pi*x)/5 + (-1)**6 * sin(6*pi*x)/6 + (-1)**7 * sin(7*pi*x)/7 + (-1)**8 * sin(8*pi*x)/8 + (-1)**9 * sin(9*pi*x)/9 )
replot 1./2 - 1/(pi)*( (-1)**1 * sin(1*pi*x)/1 + (-1)**2 * sin(2*pi*x)/2 + (-1)**3 * sin(3*pi*x)/3 + (-1)**4 * sin(4*pi*x)/4 + (-1)**5 * sin(5*pi*x)/5 + (-1)**6 * sin(6*pi*x)/6 + (-1)**7 * sin(7*pi*x)/7 + (-1)**8 * sin(8*pi*x)/8 + (-1)**9 * sin(9*pi*x)/9 + (-1)**10 * sin(10*pi*x)/10 + (-1)**11 * sin(11*pi*x)/11 + (-1)**12 * sin(12*pi*x)/12 + (-1)**13 * sin(13*pi*x)/13 + (-1)**14 * sin(14*pi*x)/14 + (-1)**15 * sin(15*pi*x)/15 + (-1)**16 * sin(16*pi*x)/16 + (-1)**17 * sin(17*pi*x)/17 + (-1)**18 * sin(18*pi*x)/18 + (-1)**19 * sin(19*pi*x)/19 + (-1)**20 * sin(20*pi*x)/20 + (-1)**21 * sin(21*pi*x)/21 + (-1)**22 * sin(22*pi*x)/22 + (-1)**23 * sin(23*pi*x)/23 + (-1)**24 * sin(24*pi*x)/24 + (-1)**25 * sin(25*pi*x)/25 + (-1)**26 * sin(26*pi*x)/26 + (-1)**27 * sin(27*pi*x)/27 + (-1)**28 * sin(28*pi*x)/28 + (-1)**29 * sin(29*pi*x)/29 + (-1)**30 * sin(30*pi*x)/30 + sin(31*pi*x)/31 )
# contemplate this for a moment and then let's simplify it by just
# plotting the last of those partial series:
reset
set samples 1000
plot 1./2 - 1/(pi)*( (-1)**1 * sin(1*pi*x)/1 + (-1)**2 * sin(2*pi*x)/2 + (-1)**3 * sin(3*pi*x)/3 + (-1)**4 * sin(4*pi*x)/4 + (-1)**5 * sin(5*pi*x)/5 + (-1)**6 * sin(6*pi*x)/6 + (-1)**7 * sin(7*pi*x)/7 + (-1)**8 * sin(8*pi*x)/8 + (-1)**9 * sin(9*pi*x)/9 + (-1)**10 * sin(10*pi*x)/10 + (-1)**11 * sin(11*pi*x)/11 + (-1)**12 * sin(12*pi*x)/12 + (-1)**13 * sin(13*pi*x)/13 + (-1)**14 * sin(14*pi*x)/14 + (-1)**15 * sin(15*pi*x)/15 + (-1)**16 * sin(16*pi*x)/16 + (-1)**17 * sin(17*pi*x)/17 + (-1)**18 * sin(18*pi*x)/18 + (-1)**19 * sin(19*pi*x)/19 + (-1)**20 * sin(20*pi*x)/20 + (-1)**21 * sin(21*pi*x)/21 + (-1)**22 * sin(22*pi*x)/22 + (-1)**23 * sin(23*pi*x)/23 + (-1)**24 * sin(24*pi*x)/24 + (-1)**25 * sin(25*pi*x)/25 + (-1)**26 * sin(26*pi*x)/26 + (-1)**27 * sin(27*pi*x)/27 + (-1)**28 * sin(28*pi*x)/28 + (-1)**29 * sin(29*pi*x)/29 + (-1)**30 * sin(30*pi*x)/30 + sin(31*pi*x)/31 )
```

Users of geogebra can simply punch in:

```
1./2 - 1/(pi)*( (-1)**1 * sin(1*pi*x)/1 + (-1)**2 * sin(2*pi*x)/2 + (-1)**3 * sin(3*pi*x)/3 + (-1)**4 * sin(4*pi*x)/4 + (-1)**5 * sin(5*pi*x)/5 + (-1)**6 * sin(6*pi*x)/6 + (-1)**7 * sin(7*pi*x)/7 + (-1)**8 * sin(8*pi*x)/8 + (-1)**9 * sin(9*pi*x)/9 + (-1)**10 * sin(10*pi*x)/10 + (-1)**11 * sin(11*pi*x)/11 + (-1)**12 * sin(12*pi*x)/12 + (-1)**13 * sin(13*pi*x)/13 + (-1)**14 * sin(14*pi*x)/14 + (-1)**15 * sin(15*pi*x)/15 + (-1)**16 * sin(16*pi*x)/16 + (-1)**17 * sin(17*pi*x)/17 + (-1)**18 * sin(18*pi*x)/18 + (-1)**19 * sin(19*pi*x)/19 + (-1)**20 * sin(20*pi*x)/20 + (-1)**21 * sin(21*pi*x)/21 + (-1)**22 * sin(22*pi*x)/22 + (-1)**23 * sin(23*pi*x)/23 + (-1)**24 * sin(24*pi*x)/24 + (-1)**25 * sin(25*pi*x)/25 + (-1)**26 * sin(26*pi*x)/26 + (-1)**27 * sin(27*pi*x)/27 + (-1)**28 * sin(28*pi*x)/28 + (-1)**29 * sin(29*pi*x)/29 + (-1)**30 * sin(30*pi*x)/30 + sin(31*pi*x)/31 )
```

As for Desmos, it has this problem where it does not handle a
pasting-in of exponents with more than one digit. So you can punch in
this expression where exponents are surrounded by `{}`

. In addition
\(\pi\) is written as `\pi`

instead of `pi`

to make the string
pasting work for Desmos.

```
1./2 - 1/(\pi)*( (-1)^{1} * sin(1*\pi*x)/1 + (-1)^{2} * sin(2*\pi*x)/2 + (-1)^{3} * sin(3*\pi*x)/3 + (-1)^{4} * sin(4*\pi*x)/4 + (-1)^{5} * sin(5*\pi*x)/5 + (-1)^{6} * sin(6*\pi*x)/6 + (-1)^{7} * sin(7*\pi*x)/7 + (-1)^{8} * sin(8*\pi*x)/8 + (-1)^{9} * sin(9*\pi*x)/9 + (-1)^{10} * sin(10*\pi*x)/10 + (-1)^{11} * sin(11*\pi*x)/11 + (-1)^{12} * sin(12*\pi*x)/12 + (-1)^{13} * sin(13*\pi*x)/13 + (-1)^{14} * sin(14*\pi*x)/14 + (-1)^{15} * sin(15*\pi*x)/15 + (-1)^{16} * sin(16*\pi*x)/16 + (-1)^{17} * sin(17*\pi*x)/17 + (-1)^{18} * sin(18*\pi*x)/18 + (-1)^{19} * sin(19*\pi*x)/19 + (-1)^{20} * sin(20*\pi*x)/20 + (-1)^{21} * sin(21*\pi*x)/21 + (-1)^{22} * sin(22*\pi*x)/22 + (-1)^{23} * sin(23*\pi*x)/23 + (-1)^{24} * sin(24*\pi*x)/24 + (-1)^{25} * sin(25*\pi*x)/25 + (-1)^{26} * sin(26*\pi*x)/26 + (-1)^{27} * sin(27*\pi*x)/27 + (-1)^{28} * sin(28*\pi*x)/28 + (-1)^{29} * sin(29*\pi*x)/29 + (-1)^{30} * sin(30*\pi*x)/30 + sin(31*\pi*x)/31 )
```

Desmos tricks for sawtooth and triangle wave:

```
Noah says:I accidentally did this last class cause I messed up the formula
Noah says:\frac{1}{2}\sum_{n=1}^{10}\frac{1}{n}\sin\left(n\pi x\right)
Noah
Noah says:\sum_{k=1}^{10}\frac{4(1-(-1)^{\left(2k-1\right)})}{\pi^{2}\left(2k-1\right)^{2}}\cos((2k-1)\pi x)
Noah says:second one is the even triangle
```

### 13.3.2. Triangle wave

Let us look at an *even* triangle wave. The Fourier series is:

```
$ gnuplot
# at the gnuplot> prompt type:
reset
set grid
set samples 1000
set xrange [-2:2]
set yrange [-1.3:1.3]
set terminal qt lw 2
plot 0 lw 1 lc 'black'
replot [-1:0] 2*x + 1 lw 2 lc 'green'
replot [1:2] 2*(x-1) - 1 lw 2 lc 'green'
replot [-2:-1] -2*x -3 lw 2 lc 'green'
replot [0:1] -2*x + 1 lw 2 lc 'green'
# now that we have a triangle wave, we try a few terms of the
# Fourier series
replot 4*(1 - (-1)**1)/(pi**2 * 1**2) * cos(1*pi*x)
replot 4*(1 - (-1)**1)/(pi**2 * 1**2) * cos(1*pi*x) + 4*(1 - (-1)**3)/(pi**2 * 3**2) * cos(3*pi*x)
replot 4*(1 - (-1)**1)/(pi**2 * 1**2) * cos(1*pi*x) + 4*(1 - (-1)**3)/(pi**2 * 3**2) * cos(3*pi*x) + 4*(1 - (-1)**5)/(pi**2 * 5**2) * cos(5*pi*x)
replot 4*(1 - (-1)**1)/(pi**2 * 1**2) * cos(1*pi*x) + 4*(1 - (-1)**3)/(pi**2 * 3**2) * cos(3*pi*x) + 4*(1 - (-1)**5)/(pi**2 * 5**2) * cos(5*pi*x) + 4*(1 - (-1)**7)/(pi**2 * 7**2) * cos(7*pi*x)
replot 4*(1 - (-1)**1)/(pi**2 * 1**2) * cos(1*pi*x) + 4*(1 - (-1)**3)/(pi**2 * 3**2) * cos(3*pi*x) + 4*(1 - (-1)**5)/(pi**2 * 5**2) * cos(5*pi*x) + 4*(1 - (-1)**7)/(pi**2 * 7**2) * cos(7*pi*x) + 4*(1 - (-1)**9)/(pi**2 * 9**2) * cos(9*pi*x)
# after contemplating each of those, put a *ton* of terms in there:
replot 4*(1 - (-1)**1)/(pi**2 * 1**2) * cos(1*pi*x) + 4*(1 - (-1)**3)/(pi**2 * 3**2) * cos(3*pi*x) + 4*(1 - (-1)**5)/(pi**2 * 5**2) * cos(5*pi*x) + 4*(1 - (-1)**7)/(pi**2 * 7**2) * cos(7*pi*x) + 4*(1 - (-1)**9)/(pi**2 * 9**2) * cos(9*pi*x) + 4*(1 - (-1)**11)/(pi**2 * 11**2) * cos(11*pi*x) + 4*(1 - (-1)**13)/(pi**2 * 13**2) * cos(13*pi*x) + 4*(1 - (-1)**15)/(pi**2 * 15**2) * cos(15*pi*x) + 4*(1 - (-1)**17)/(pi**2 * 17**2) * cos(17*pi*x) + 4*(1 - (-1)**19)/(pi**2 * 19**2) * cos(19*pi*x) + 4*(1 - (-1)**21)/(pi**2 * 21**2) * cos(21*pi*x) + 4*(1 - (-1)**23)/(pi**2 * 23**2) * cos(23*pi*x) + 4*(1 - (-1)**25)/(pi**2 * 25**2) * cos(25*pi*x) + 4*(1 - (-1)**27)/(pi**2 * 27**2) * cos(27*pi*x) + 4*(1 - (-1)**29)/(pi**2 * 29**2) * cos(29*pi*x) + 4*(1 - (-1)**31)/(pi**2 * 31**2) * cos(31*pi*x) + 4*(1 - (-1)**33)/(pi**2 * 33**2) * cos(33*pi*x) + 4*(1 - (-1)**35)/(pi**2 * 35**2) * cos(35*pi*x) + 4*(1 - (-1)**37)/(pi**2 * 37**2) * cos(37*pi*x) + 4*(1 - (-1)**39)/(pi**2 * 39**2) * cos(39*pi*x) + 4*(1 - (-1)**41)/(pi**2 * 41**2) * cos(41*pi*x) + 4*(1 - (-1)**43)/(pi**2 * 43**2) * cos(43*pi*x) + 4*(1 - (-1)**45)/(pi**2 * 45**2) * cos(45*pi*x) + 4*(1 - (-1)**47)/(pi**2 * 47**2) * cos(47*pi*x) + 4*(1 - (-1)**49)/(pi**2 * 49**2) * cos(49*pi*x) + 4*(1 - (-1)**51)/(pi**2 * 51**2) * cos(51*pi*x) + 4*(1 - (-1)**53)/(pi**2 * 53**2) * cos(53*pi*x) + 4*(1 - (-1)**55)/(pi**2 * 55**2) * cos(55*pi*x) + 4*(1 - (-1)**57)/(pi**2 * 57**2) * cos(57*pi*x) + 4*(1 - (-1)**59)/(pi**2 * 59**2) * cos(59*pi*x)
```

### 13.3.3. The Gaussian

The Fourier transform of a gaussian is fascinating from more than one point of view: the mathematical tricks used to calculate it are delightful and the applications are metaphorically rich.

#### 13.3.3.1. Visuals on the Gaussian

Let us start by exploring the Gaussian function. We will use the form
that gives the *probability density function* for a normally
distributed random variable \(x\) with mean \(\mu\) and
standard deviation \(\sigma\):

If our mean is zero and standard deviation is one then we get:

The Gaussian function comes up in so many places that it pays to become very comfortable with its behavior. Let us start with some plots to get a feeling for what it looks like:

```
reset
set grid
set terminal qt lw 3
set samples 1000
set xrange [-5:5]
set yrange [0:1]
plot (1.0/(1*sqrt(2*pi))) * exp(- (1./2) * (x - 0)**2 / 1**2) title 'mean 0 sigma 1'
replot (1.0/(2*sqrt(2*pi))) * exp(- (1./2) * (x - 0)**2 / 2**2) title 'mean 0 sigma 2'
replot (1.0/(0.5*sqrt(2*pi))) * exp(- (1./2) * (x - 0)**2 / 0.5**2) title 'mean 0 sigma 0.5'
replot (1.0/(0.4*sqrt(2*pi))) * exp(- (1./2) * (x - 1.3)**2 / 0.4**2) title 'mean 1.3 sigma 0.4'
```

Looking closely at the formula for the Gaussian you will notice that \(\sigma\) appears in the denominator of the exponential, as well as the denominator of the constant scale factor \(\frac{1}{\sigma\sqrt{2\pi}}\). This means that:

- \(\sigma\) is big
The curve will spread out and not be tall. This matches our intuition: if there is a big uncertainty then our random variable will often be quite different from the mean value.

- \(\sigma\) is small
The curve will be skinny and taller. The intuition here is that random values are likely to be close to the mean, and not

*deviate*from it much.

#### 13.3.3.2. Normalization and total probability

If we think of the Gaussian as a probability distribution function then it is important that the probabilities of a random variable being in all its possible states add up to 1. This means that we insist that this integral:

equal 1 for all values of \(\mu\) and \(\sigma\).

This is guaranteed by the *Gaussian integral* formula:

If we substitute \(x \rightarrow ax\) then we get:

For the Gaussian our \(a^2\) is \(1 / (2\sigma^2)\), so we get:

which implies that

This means that the Gaussian has been crafted carefully so that, no matter what \(\sigma\) is, we always end up with a total probability of 1.

The factor of \(\frac{1}{\sigma\sqrt{2\pi}}\) which multplies our
exponential function is called a *normalization constant*.
Normalization constants are widely used term in mathematics: we
apply a constant scale to a function so that it satisfies some
necessary condition.

Note that I have not mentioned the mean in discussing these integrals. The reason is that the gaussian goes from \(-\infty\) to \(\infty\), so a shift of the hump to the right or to the left does not affect the integral at all.

#### 13.3.3.3. Fourier-transforming the Gaussian

The first thing I have to do here is be a bit apologetic: the Gaussian
is *not* a periodic function. Nor is it defined on a finite domain
(which can be treated as a periodic function. Rather, it unfolds on
the entire real line with no repetitions.

This means that we cannot apply the mechanisms we saw earlier to calculate the Fourier transform of the Gaussian: this will be a different beast.

But the beast does exist: you can take an arbitrary function and write
it down as a *Fourier integral* (not a sum; an integral!) which
integrates sin and cos functions with a continuum of frequencies, and
we will get our function back.

The Fourier transform of this more general function \(f(t)\) will
not be a collection of discrete *coefficients* \({a_k}, {b_k}\)
(remember: those \(k\) indices were multiples of a fundamental
frequency).

Instead it will involve a pair of *functions* \(A(\omega),
B(\omega)\) where \(\omega\) is the angular frequency \(\omega
= 2 \pi f\):

You then calculate \(A(\omega)\) and \(B(\omega)\) quite analogously to the discrete Fourier coefficients:

Let us apply this to the Gaussian. We will center it at a mean of zero then we have an even function. This means that \(B(\omega) = 0\), so we have: transform:

the application of some integration techniques gives us:

The important thing to take home here is that *the Fourier transform
of a Gaussian is another Gaussian!* And this other Gaussian has
\(2/\sigma\) in all the places where the original one has
\(\sigma\). If you let \(\delta = 2/\sigma\) then you get:

so we can say that:

The Fourier transform of a Gaussian is a gaussian with the reciprocal standard deviation.

Caution

The equations above for \(f(t)\) and \(A(\omega)\) are subject to normalization. Researchers in different fields use a different normalization convention, with factors of 2 and \(\pi\) (possibly with square roots) that can be applied in one formula or in the other, as long as it is consistent. The lack of a single convention is frustrating, and ends up being costly, as collaborators have to be very careful. At this time I have not yet double-checked that I am using a uniform normalization convention in the equations above. Still, the fundamental conclusion holds: the transform of a Gaussian is a Gaussian with reciprocal standard deviation.

#### 13.3.3.4. Visualization

Here is an animated visualization of the Gaussian and its transform for values of sigma that go from 10 to 1/10.

```
set grid
set samples 1000
do for [ii=1:100] {
sigma = 10.0 / ii
plot 1./(sigma*sqrt(2*pi)) * exp(-x**2 / (2*sigma**2))
replot exp(-x**2 * sigma**2 / 4)
pause(1.3)
}
```

We will show more visualizations when we look at uses of the Fourier transform on discrete data.

#### 13.3.3.5. Thoughts about the transform of the Gaussian

[this section is not yet elaborated]

The general idea that if you have a narrowly localized signal in time, its Fourier transform will be spread out in frequency. Conversely, if you have a spread out signal, then its Fourier transform will be concentrated in a narrow band of frequencies.