Fourier analysis is the study of how a function can be represented as a sum of waves. Take a look at the animation playing at the side, a shape is being drawn using a chain of rotating circles of different sizes. You can even try drawing your own shape, it's interactive! This article will cover in detail how this animation works, and what math is behind it. The concepts that we'll discover are widely used in electronics, acoustics and communications. Operators such as the Fourier Transform are constantly used in the real world, without these discoveries the world would not be the same. Much software relies in Fourier Analysis, such as for instance Shazam, the famous service for identifying songs. Any audio spectrum visualized processes the signal using Fourier Transform, these are just a few of the many application of this analysis.
This article requires a general good understanding about math. If you want to fully understand every bit of math, there are a few requirements:
Even if you don't understand the math involved, the main ideas about Fourier analysis will be explained visually and occasionally with interactive examples.
If you understand the following expressions and notations you should have a solid background
Mathematician and physicist Joseph Fourier determined that a function can be represented as a series of sines and cosines of different frequencies and different amplitudes.
Fourier is notoriously known for having developed Fourier series and the Fourier transform, which are the main focus of this article.
We will also cover the FFT algorithm, which is a fast implementation of the Fourier transform widely used among software.
Jean Baptiste Joseph Fourier (Auxerre, 1768 - Paris, 1830)
The Fourier series is the representation of a periodic function with a summation of sine and cosine waves of discrete frequencies. Each wave is weighted according to "how important" it is to represent the original function.
Fourier Series are often represented in two ways: trigonometric and exponential. They both work in the same way, but the exponential one is also defined on the complex plane and as we'll see, has a nicer, more elegant form.
The Fourier transform is an operation that transforms a signal from time-domain to a continuous frequency-domain. The function can be a generic, not necessarily period function \(f(x)\). The output of the Fourier transform \(\mathcal{F}\) is a complex-valued function whose absolute value represents the magnitude of each frequency.
$$ \mathcal{F}\{f(t)\}=\hat{f}(\xi) $$A function \(f(t)\) is periodic if there is a positive number \(T\) (the period of \(f\)) such that
$$ f(t+nT)=f(t)\quad \forall t \in D_f, n \in \mathbb{Z} $$We can represent a periodic function using a sum of sines and cosines, for each discrete frequency we have a wave with its own weight (its amplitude).
$$ f(t)=C+\sum_{n=1}^{\infty} a_n \cos\left(\frac{2\pi nt}{T}\right) + b_n \sin\left(\frac{2\pi nt}{T}\right) $$We're computing a sum from \(n=1\) to infinity, where \(n\) represents each discrete frequency. The term \(\frac{2\pi n}{T}\) controls the frequency based on \(n\). A normal sine of cosine wave oscillates every \(2\pi\), with this modification e.g. \(n=5\) means that the function will oscillate \(5\) times within that span. The terms \(a_b\) and \(b_n\) control how much that particular frequency is important. You might notice that with this we can only represent functions "laying" on the \(x\)-axis. To resolve this problem, we add a generic constant term \(C\) to shift the function up/down.
>> Note: We only need the discrete frequencies (1 Hz, 2 Hz, ...) to represent the function, although the Fourier transform gives us a continuous frequency analysis.
Sum of waves
Great! One small problem, what are \(a_n\), \(b_n\) and \(C\)? Here we manipulate the Fourier series definition to find these values given \(f(t)\) on a period \([t_0;t_0+T]\).
>> Note: To simplify the equations i'm going to define \(w_k=\frac{2\pi k}{T}\)
Starting from the \(C\) term, we take the integral over the generic period \([t_0;t_0+T]\) on both sides
$$ \integral{t_0}{t_0+T}{f(t)}{t} = \integral{t_0}{t_0+T}{h}{t} + \sum_{n=1}^{\infty} \left[ a_n \integral{t_0}{t_0+T}{\cos(w_n t)}{t} + b_n \integral{t_0}{t_0+T}{\sin(w_n t)}{t} \right] $$If you think about it, the integral over a full period of a function such as \(\sin(x)\) or \(\cos(x)\) is 0. If we consider \(\sin(w_n x)\) or \(\cos(w_n x)\) the function will make more full cycles in the span of the interval \(T\), all of which yield an area of 0.
$$ \begin{align*} \integral{t_0}{t_0+T}{f(t)}{t} &= \integral{t_0}{t_0+T}{C}{t} + \sum_{n=1}^{\infty} a_n\cdot 0 + b_n\cdot 0\\ &= \integral{t_0}{t_0+T}{C}{t} \\ &= C \integral{t_0}{t_0+T}{}{t} \\ &= C {\left[x\right]}^{t_0+T}_{t_0} \\ &= C(t_0+T-t_0) \\ &= CT \end{align*} $$ concluding that $$ C = \frac{1}{T} \integral{t_0}{t_0+T}{f(t)}{t} $$Again, we take the integral over the period \([t_0;t_0+T]\) on both sides, but first we multiply everything by \(\cos(w_k t)\)
$$ \begin{align*} \integral{t_0}{t_0+T}{f(t)\cos(w_k t)}{t} &= h \integral{t_0}{t_0+T}{\cos(w_k t)}{t} + \sum_{n=1}^{\infty} \left[ a_n \integral{t_0}{t_0+T}{\cos(w_n t)\cos(w_k t)}{t} + b_n \integral{t_0}{t_0+T}{\sin(w_n t)\cos(w_k t)}{t} \right] \end{align*} $$By using orthogonality relationships or by literally evaluating the above integrals, we get the following
$$ \begin{align*} \integral{t_0}{t_0+T}{f(t)\cos(w_k t)}{t}&= \integral{t_0}{t_0+T}{a_k\cos^2(w_k t)}{t} \\ &= a_k \left(\frac{T}{2}\right) \end{align*} $$concluding that
$$ a_k=\frac{2}{T}\integral{t_0}{t_0+T}{f(t)\cos(w_k t)}{t} $$To find \(b_n\), we do the exact same thing, but instead of multiplying by \(\cos(w_n t)\), we multiply by \(\sin(w_n t)\), ending up with
$$ b_k=\frac{2}{T}\integral{t_0}{t_0+T}{f(t)\sin(w_k t)}{t} $$>> Note: \(C=\frac{a_0}{2}\)
Under appropriate conditions, we can represent a periodic function \(f(t)\) on an interval \(t_0;t_0+T\) by
$$ \begin{align*} f(t)&=\frac{a_0}{2}+\sum_{n=1}^{\infty} a_n \cos(w_n t) + b_n \sin(w_n t) \\ a_k &= \frac{2}{T} \integral{t_0}{t_0+T}{f(t)\cos(w_k t)}{t} \\ b_k &= \frac{2}{T} \integral{t_0}{t_0+T}{f(t)\sin(w_k t)}{t} \end{align*} $$where
$$ w_k=\frac{2\pi k}{T} $$To start, we will plot our signal in the complex plane. We will plot it while making it rotate around the origin \(0+0i\). Recall that Euler's Formula tells us that the function \(e^{it}\) is a rotation around the origin, also called the unit circle. If we multiply our signal by this circle, it will follow its path, achieving a polar plot-like graph: \(f(t)e^{it}\). Now, the rotational function makes a full cycle every \(2\pi\), we can plot our signal at a different speed (different frequencies). To do so, we will multiply the time of the rotational function by \(2\pi \xi\) where \(\xi\) is the frequency. We add the \(2\pi\) term so that if your frequency is \(1\), we will have a rotation each second rather then every \(2\pi\) seconds \((\sim 6.28\,s)\). Furthermore, we want the unit circle to rotate clockwise instead of counter-clockwise. We will just add the negative sign to the time \(t\) argument of the rotational function. Our final function for now is \(f(t)e^{-2\pi it\xi}\). You can see the animation next to this paragraph, you can vary the frequency or even draw your own signal.
The next step is to create a new function. Instead of plotting our signal as time increases with a given frequency \(\xi\), we plot all of the signal at once but the frequency changes over time. Here we are moving from a time-dependant function to a frequency-dependant function. The argument is no longer the time \(f(t)\) but rather the frequency with which we want to plot our signal around the origin \(f(\xi)\). Now, what is the center of mass? The center of mass is basically the average point of the function, which is a complex number. You might notice that the center of mass (the blue dot) in the animation changes as the frequency changes. To find its value for a certain frequency, we can sample a bunch of points from the function and then divide it by the number of samples. This approach works when we are dealing with discrete functions, such as the animation (the signal you draw is a set of points), but we would need an infinite number of samples when the signal is continuous. An infinite amount of precision is achievable using an integral. To find the center of mass when the frequency is \(\xi\) we integrate our complex function over a certain period of time, and then divide it by the time length.
$$ \frac{1}{t_1-t_0} \integral{t_0}{t_1}{f(t)e^{-2\pi i t\xi}}{t} $$The next animation is distance of the center of mass from the origin as the frequency changes. The key concept is that when the frequency matches the period, the center of mass is unusually further from the origin. Go back to the last section and look at the blue dot moving far away from the origin when the frequency is the right one (reset the sine or cosine wave so that it is clearer). When the frequency is a component of the function, this distance peaks. If the function is composed of multiple frequencies, the same distance will peak multiple times, and more it peaks, the more the frequency is present in the original function.
The function that represents the center of mass as the frequency changes is called Fourier transform. Actually, not quite, the Fourier transform is defined with the integral over \(\mathbb{R}\) and is not divided by the time span. The absolute value (distance from the origin) of this function is the amount of all continuous frequencies present in the original signal.
$$ \hat{f}(\xi)=\integral{-\infty}{\infty}{f(t)\,e^{-2\pi it\xi}}{t} $$Let's look at a simple example. We are going to derive the Fourier series of a function \(f(x)\) defined as such:
$$ f(x)= \begin{cases} -1\quad \text{if } 0 \lt x \lt \pi \\ +1\quad \text{if } -\pi \lt x \lt 0 \\ \end{cases} $$The period of this function is \(T=2\pi\). We can already simplify the \(\frac{2\pi}{T}\) term, leaving us with
$$ f(x)=\frac{a_0}{2} + \sum_{n=1}^{\infty} a_n \cos(nx) + b_n \sin(nx) $$First, we need to find \(a_n\). Simplifying \(\frac{2\pi}{T}\) and \(\frac{T}{2}\) we get
$$ a_n=\integral{-\pi}{\pi}{f(x)\cos(nx)}{x} $$Looking at the graph we notice that we can split the integral into two parts at \(x=0\). On the left part, the function is \(-\cos(nx)\), while on the right part the function is \(\cos(nx)\).
$$ \begin{align*} a_n &= \frac{1}{\pi} \integral{-\pi}{0}{-\cos(nx)}{x} + \frac{1}{\pi} \integral{0}{\pi}{\cos(nx)}{x} \\ &= -\frac{1}{\pi} \integral{-\pi}{0}{\cos(nx)}{x} + \frac{1}{\pi} \integral{0}{\pi}{\cos(nx)}{x} \\ &= -\frac{1}{\pi} {\left[\frac{\sin(xn)}{n}\right]}_{-\pi}^{0} + \frac{1}{\pi} {\left[\frac{\sin(xn)}{n}\right]}_{0}^{-\pi} \\ &= -\frac{1}{\pi} \left[\frac{\sin(\pi n)}{n}\right] + \frac{1}{\pi} \left[\frac{\sin(\pi n)}{n}\right] \\ &= \left(\frac{1}{\pi}-\frac{1}{\pi}\right) \left[\frac{\sin(\pi n)}{n}\right] \\ &= 0 \end{align*} $$\(a_n\) is always going to be 0. We can remove the \(a_n \cos(nx)\) and \(\frac{a_0}{2}\) terms from the series.
Now the same thing for \(b_n\)
$$ b_n=\integral{-\pi}{\pi}{f(x)\sin(nx)}{x} $$Again, we split the integral into two parts
$$ \begin{align*} b_n &= -\frac{1}{\pi} \integral{-\pi}{0}{\sin(nx)}{x} + \frac{1}{\pi} \integral{0}{\pi}{\sin(nx)}{x} \\ &= -\frac{1}{\pi} {\left[\frac{-\cos(xn)}{n}\right]}_{-\pi}^{0} + \frac{1}{\pi} {\left[\frac{-\cos(xn)}{n}\right]}_{0}^{-\pi} \\ &= -\frac{1}{\pi} \left[-\frac{1}{n}+\frac{\cos(\pi n)}{n}\right] + \frac{1}{\pi} \left[\frac{-\cos(\pi n)}{n}+\frac{1}{n}\right] \\ &= -\frac{1}{\pi} \left[\frac{\cos(\pi n)-1}{n}\right] + \frac{1}{\pi} \left[\frac{1-\cos(\pi n)}{n}\right] \\ &= \frac{2}{\pi} \cdot \frac{1-\cos(\pi n)}{n} \\ &= \frac{2-2\cos(\pi n)}{\pi n} \end{align*} $$Given \(b_n\) our series is now complete!
$$ f(x)=\sum_{n=1}^{\infty} \frac{2-2\cos(\pi n)}{\pi n} \cdot \sin(nx) $$We won't simplify this further, therefore this is our final result.
The effort pays off when we graph this function, as more terms are added, the function looks more and more like the original square wave. You can drag the slider by opening the left panel.
>> Note: If \(f(x)\) is even, the coefficient \(b_n\) will always be equal to zero. If \(f(x)\) is odd, \(a_n\) will always be equal to zero.
The exponential Fourier Series is the same thing but extended to the complex plane, we can use Euler's identify to manipulate the real Fourier series and get the following expression
$$ f(t)=\frac{a_0}{2}+\sum_{n=1}^{\infty} \frac{1}{2}(a_n-ib_n)e^{iw_n t}+ \frac{1}{2}(a_n+ib_n)e^{-iw_n t} $$Now, let \(n\) be also negative. With an appropriate coefficient \(c_n\) we can arrange the series as
The coefficient \(c_n\) can be computed as
$$ c_n = \frac{1}{T} \integral{t_0}{t_0+T}{f(t)e^{-itw_n}}{t} $$The Fast Fourier Transform is a computer algorithm to compute the Discrete Fourier Transform. The fastest and most used implementation of this algorithm is FFTw, a C subroutine library written at MIT.
When computing the DFT on a (discrete) signal \(f(t)\), we take the average of all the points of \(f(t)e^{-2\pi i\xi}\). This process is repeated for each frequency \(\xi\). Computing all of these values is a \(O(N^2)\) operation, but with the FFT algorithm we can decrease the number of operations to \(O(N\log(N))\). Many FFT algorithms depend on that fact that \(e^{-2\pi i/N}\) is a root of unity.
There are plenty of FFT algorithms, here's a few: Cooley–Tukey FFT, Prime-factor FFT, Bruun's FFT, Rader's FFT, Chirp Z-transform, Hexagonal fast Fourier transform
There is also another version called SFT (Sparse Fourier Transform), which is a DFT for handling big data signals, and is mainly used in GPS synchronization.
So how does the cool epicycles animation work? First of all we need apply the Fourier transform operation, however the drawing is just a set of points, it's a discrete function rather than a continuous one, this means that we'll need to use the Discrete Fourier Transform operator. Each circle represents a discrete frequency, and each center revolves around the previous circle's circumference. The radius of the circle is the magnitude of the current frequency, which is the absolute value of the Fourier transform at that frequency. The revolution is based on the time passed and the phase of the Fourier transoform.