Fourier Analysis (Introduction)

Discrete Fourier Transformation

Emre S. Tasci emre.tasci@hacettepe.edu.tr
01/09/2021

$\newcommand{\unit}[1]{\,\text{#1}}$

Nature hosts a collection of harmonious processes and we do our best to imitate it. Take the music, for example, or paintings... or the perfect ordering of crystals. At our hands, we have the almost infinite possibilities though we can't say we harness most of it. When a friend asks about our favorite number, there is a very slim chance that our answer will include a number greater than 100 or for that matter a negative or even an irritational or a complex one ("3", I hear you saying 8). We confine ourselves to within a sigma of a "normal" values -- the dominancy of Gaussian distribution in almost everything we touch is a proof of that. But this is neither a boring, nor a fruitless trait, actually there is a good side to it as well, this way we can analyze things much more easily and one way to do it is the Fourier Analysis: the decomposition of complex looking signals, functions (,information) into simple sinusoids. In this treatise we'll see how it works and is applied to.

This work is limited to discrete Fourier transformation. The good news is almost all of our data comes as discrete points, from a number of measurements (or calculations).

We will start by defining the characteristics of a sinusoids, that of a wave to put it simply.

Period ($T$), Frequency ($f$) and Angular Frequency ($\omega$)

Consider a signal such as the sound wave or the shape of a rope while playing skipping rope. We immediately match their corresponding units: for the sound wave, it is the air pressure (manifested for example as the displacement of the eardrum from its equilibrium position) at a given time; while for the rope, it is the displacement of a segment of the rope at a given position. But in essance, our function representing the signal may look like this:

Looking at the formula, $\sin(2\pi \frac{t}{3})$ we observe two additional coefficients next our parameter $t$: $2\pi$ and 3. "3" is the period of our signal ($T$). It means that, it takes 3 seconds to complete 1 sine wave. We have taken our observation interval from t = 0 to t = 12 s, therefore, we observe 12 / 3 = 4 full waves.

In this example, we were 'lucky' to take our observation time to 12 s, as it was an integer multiple of our period of 3 s, so we were able to directly say "4 full waves in 12 seconds". But what about the question of how many waves had passed by the 5th second? We can resort to proportionality: if in $T$ seconds 1 wave passes, then how many waves pass in 1 second? The answer is obviously: $\tfrac{1}{T}$ waves, which is defined as the frequency ($f$): the number of waves passing per second. Its unit is $\tfrac{1}{\unit{s}} = \unit{s$^{-1}$} = \unit{Hz}$ if we are measuring the change of a signal with respect to time (like a sound wave's evolution wrt time), and $\unit{m$^{-1}$}$ if we are measuring the change wrt length (like that of a rope being oscillated).

Checking the definitions, we concur that the period and the frequency are related as: $$f = \frac{1}{T} \leftrightarrow T = \frac{1}{f}$$

In general, we can write a sine or a cosine wave in the following form: $$f(t) = A\sin\left(2\pi\frac{t}{T} + \phi \right)$$

Here $A$ is the amplitude, that enables to get values outside of $[-1,1]$ domain of the sine function and $\phi$ is the phase that shifts the signal left and right (remember that cosine is just sine shifted $\tfrac{\pi}{2}$ towards left. Keeping in mind that we can always write cosine as sine (and vice versa) by introducing a phase, we'll drop the phase and from now on will be working with cosines.

$$f(t) = A\cos\left(2\pi\frac{t}{T}\right) = A\cos\left(\frac{2\pi}{T}t\right) =A\cos\left(2\pi f t\right)$$

When we defined/derived frequency, we asked about the number of waves that passed in 1 second. But we can also think of the repeating motion as a rotation around a circle: at t = 0, let's say we start moving from $(1,0)$ point on the unit circle with a constant speed, then at t = T/4 we will be passing from $(0,1)$ point, at t = T/2 by the $(-1,0)$, t = 3T/4 at $(0,-1)$ and finally at t = T back at $(1,0)$. As illustrated in this image:

HarmonicMotion
Harmonic Motion (from Wikipedia)

The point and beauty is, the radius of the circle actually does not matter: instead of the unit circle, it could have been a much bigger circle but given that we complete our rotation around it in T seconds, we would still be on the "top" ($\theta=\tfrac{\pi}{2}$) at t = T/4, on the "left" ($\theta=\pi$) at t = T/2, at the "bottom" ($\theta=\tfrac{3\pi}{2}$) at t = 3T/4, and back at our starting point ($\theta=2\pi$) at t = T, so it may prove advantageous to be asking the question of "how many radians are covered in 1 second?" If it takes T seconds to complete $2\pi$, then the answer is simply $\tfrac{2\pi}{T}$ radians per second (or substituting $f$ for $T$: $2\pi f$). This quantity is defined as the angular frequency ($\omega$).

$$f(t) = A\cos\left(\frac{2\pi}{T}t\right) =A\cos\left(2\pi f t\right)=A\cos\left(\omega t\right)$$

and as can be easily seend, the angle covered in $t$ seconds is just $\omega t$ (radians).

Sampling

When we are making a measurement, we are never measuring the whole, continous signal - we can make as many observations as we want but it is always the points from the "invisible" signal that we are observing. The more observations we make, the more reliable our graph will be with respect to the actual signal but nevertheless, always discrete. This is even present in the graphs we plot in Python.

Consider plotting a complete sine wave using 10 points (for practicality, let's assume they are uniformly distributed):

This is how we would perceive the signal if we had taken only 10 measurements. Increasing the number of observations, of course, takes us closer to the "ideal" wave that's being broadcasted. Here is another take with $N = 50$ samples:

and here's the one with $N=500$ samples:

so, the number of data that we use to identify the characteristics of the material/signal/"thing" is crucial.

Let's redefine another signal in the shape of a cosine wave, values alternating between $\pm3.2$ and taking $T=1/3\unit{s}$ to complete 1 wave. Hence:

$$x(t) = 3.2\cos(2\pi (3) t)$$

The $3$ factor in the equation above is the frequency, i.e., $f=\frac{1}{T}=\frac{1}{1/3} = 3$.

Observing this wave for 4 seconds by making a total of 250 measurements over this duration, we'll see something like this:

Let's talk about the sampling: we didn't explicitly say but it is 250 measurements uniformly distributed over the 4 seconds (as defined by the linspace function), so they will be spaced by $\tfrac{4}{250}\unit{s} = 0.016\unit{s}$ apart from each other.

The above plot is actually not a continous line but these discrete 250 points connected by lines between them:

Let's zoom in, say between t = 1.75 s to t = 2.00 s

So this is what sampling is: number of points taken (observations/measurements made) in a given interval. These measurements can be done uniformly like we did above or in a random way as well.

Let's do a random sampling of 250 points:

As the tr points, i.e. the instances at when we made our measurements, are not ordered, connecting them sequentially with lines creates such a mess. We'll have more luck in just plotting the data without connecting them.

Looks like our old cosine function. To see it's indeed that, let's sort the tr points and their corresponding values. We can manage this task two ways:

  1. Sorting tr and then recalculating xr, or,
  2. Sorting tr and re-arranging xr accordingly

The 1st approach is much more practical for this case (and also, we should have done it so in the first time!) as the calculation of 250 points of a very simple function would take no time.

However, in real cases, the evaluation of a function can be costly and thus, re-arranging accordingly option would be much more beneficial. For this purpose, we'll follow the 2nd path and will employ np.argsort() function which returns the ordered indices of the sorted array.

At some points, discrepancies can be seen but not that bad after all.

What would we end up with if instead of 250 samples, we had made 50 observations?

(such a disaster! 8)

so, the moral of story is: Sampling matters!

_even though at the moment random sampling does not seem to be a good idea, we'll show that it actually is (on the long run) - good enough to even beat the Nyquist rate itself! (more on that in the future, and certainly with a couple of accompanying conditions 8)_

We will return to the importance of sampling once again when we get a little bit into the Fourier transformations.

Discrete Fourier Transformation (DFT*)

We can define the signal we used above by either using the -ideally an infinite number of- samples, or by stating its amplitude & frequency. In our case, the 1st description method involves 250 values, hence the function is represented by a 250-dimensional vector:

The second approach is much more data-space saving with just the following 3 information:

The two approaches -in principle- yield the same signal. For some cases the first representation takes less data-space (imagine a signal that has a value of 3.2 at t = 0 point and 0 everywhere else), and for the others (like our case) the second one.

Fourier transformation is a way to go from one representation to the other one. As our data are discrete, we will be working with the Discrete Fourier Transform, defined as:

$$\begin{align*} X_k &= \sum_{n=0}^{N-1} x_n \cdot e^{-i\frac { 2\pi}{N}kn}\\ &= \sum_{n=0}^{N-1} x_n \cdot \left[\cos\left(\frac{2 \pi}{N}kn\right) - i \cdot \sin\left(\frac{2 \pi}{N}kn\right)\right] \end{align*}$$

To go from $\vec X$ to $\vec x$, we can employ its reciprocal formulation:

$$\begin{align*}x_n &= \frac{1}{N} \sum_{k=0}^{N-1} X_k \cdot e^{\,i\frac{ 2 \pi}{N}kn}\\ &= \frac{1}{N}\sum_{k=0}^{N-1} X_k \cdot \left[\cos\left(\frac{2 \pi}{N}kn\right) + i \cdot \sin\left(\frac{2 \pi}{N}kn\right)\right] \end{align*} $$

Among many fantastic properties of Fourier transformations, two are of critical importance to us:

  1. Fourier transformation is a 1:1 transformation. As a consequence: if $\vec x$ have $N$ elements, then its transformation $\vec X$ will also be of $N$ elements.
  2. Fourier transformation is a linear transformation, that is, if $\vec X$ and $\vec Y$ are the Fourier transformations of $\vec x$ and $\vec y$, respectively, then the Fourier transformation of $a \vec x + b\vec y$ will be $a\vec X + b\vec Y$. ($a,b\in \mathbb{C}$).

* [Side note: there is a very effective and widely used method in physics called the _Density Functional Theory_ which is used to solve the Schrödinger equation, also abbreviated as "DFT", so when you see the letters DFT do not immediately assume that it stands for the Discrete Fourier Transformation!]

"Manual" Application

Before we proceed any further, let's -for a moment- forget about the formal definition of Fourier transformation and instead, focus on writing our signal in terms of addition of functions in the form of:

$$ \sum X_k \left[\cos(2\pi f_k t) + i \sin(2\pi f_k t)\right]$$

since $\{X_k\}$ coefficients are complex, let's write it as $(a_k+ib_k)$:

$$ (a_k+ib_k) \left[\cos(2\pi f_k t) + i \sin(2\pi f_k t)\right]$$

distributing the coefficient inside the brackets yields:

$$\begin{align*}&\sum (a_k+ib_k) \left[\cos(2\pi f_k t) +i \sin(2\pi f_k t)\right]=\\ &=\sum\left\{ a_k\cos(2\pi f_k t) + ia_k\sin(2\pi f_k t)+ib_k\cos(2\pi f_k t)-b_k\sin(2\pi f_k t)\right\}\\ &=\sum\left\{\left[a_k\cos(2\pi f_k t) -b_k\sin(2\pi f_k t)\right]+i\left[b_k\cos(2\pi f_k t)+a_k\sin(2\pi f_k t)\right]\right\}\\ &=\sum\left[a_k\cos(2\pi f_k t) -b_k\sin(2\pi f_k t)\right]+i\sum\left[b_k\cos(2\pi f_k t)+a_k\sin(2\pi f_k t)\right] \end{align*}$$

A very difficult question...

Now comes a difficult question: to represent our signal that was $\boxed{x(t) = 3.2\cos(2\pi (3) t)}$, how should we assign values to $\{a_n\}$ and $\{b_n\}$ coefficients while keeping $n$ as small as possible (in other words, to include as less contributing functions as possible)?

( ... please wait while the reader ponders ... )

The first thing that may come to your mind could be to take $a_0 = 3.2$, $b_0 = 0$ and $f_0 = 3$, but immediately after you'll realize that it won't work because the formula forces us to take the $b_n\sin(\dots)$ component alongside as well, so we'd end up with:

$$\mathscr{F}_0 \equiv 3.2\cos(2\pi 3 t) + i3.2\sin(2\pi 3 t)$$

and even worse is we have a complex result in our hands whereas our original signal was purely real. We need to get rid of the imaginary part containing the sine: if we can do it, we have already pinned down the cosine, so we'd be be done!

The easiest way to get rid of it would be, to define a second contributing function, this time with $a_1 = 3.2$, $b_1 = 0$ and $f_1 = -3$ (yes, I'm aware I'm defining a "negative"(!) frequency! - more on that one in a minute). Thus we'll have from this:

$$\mathscr{F}_1 \equiv 3.2\cos(2\pi (-3) t) + i3.2\sin(2\pi (-3) t)$$

Applying the trigonometric identities $\cos(-\theta) = \cos(\theta)$ and $\sin(-\theta) = -\sin(\theta)$ we get:

$$\mathscr{F}_1 \equiv 3.2\cos(2\pi (-3) t) + i3.2\sin(2\pi (-3) t) = 3.2\cos(2\pi 3 t) - i3.2\sin(2\pi 3 t) $$

Almost there! So we used to contributing factors from the Fourier formula, and if sum them up now, let's see what happens:

$$\begin{align*} \mathscr{F}_0 + \mathscr{F}_1 &= \underbrace{ 3.2\cos(2\pi 3 t) + i3.2\sin(2\pi 3 t)}_{\mathscr{F}_0}+\underbrace{3.2\cos(2\pi 3 t) - i3.2\sin(2\pi 3 t)}_{\mathscr{F}_1} \\ &=6.4\cos(2\pi 3 t)\end{align*}$$

So close! We're off only by the amplitude being double the original one (due to one coming from $\mathscr{F}_0 $, the other from $\mathscr{F}_1$ , but surely this is not a problem (it is actually the reason we have the $\frac{1}{N}$ factor standing in front of the summation symbol for the inverse Fourier transformation). If we take $a_0 = a_1=\frac{3.2}{2}=1.6$, $b_0 = b_1 =0$ and $f_0 =3$, $f_1 = -3$, we regain our original signal:

$$\begin{gather*}\mathscr{F}_0 \equiv 1.6\left[\cos(2\pi 3 t) + i\sin(2\pi 3 t)\right]\\ \mathscr{F}_1 \equiv 1.6\left[\cos(2\pi (-3) t) + i\sin(2\pi (-3) t)\right]\\ \Rightarrow \mathscr{F}_0 + \mathscr{F}_1 = 3.2\cos(2\pi 3 t) = x(t)\,\checkmark \end{gather*}$$

YES!!! Congratulations! You have managed to evalaute all the complex looking summations with trigonometric multiplications by hand!

HW: What if, instead of the cosine signal, we had a sine signal like $x(t) = 3.2\sin(2\pi3t)$? What would be the coefficients then? Can you derive them without evaluating the summations and all?

Evaluating the summations

Let's have the computer a go at the tedious job of evaluating the summations to calculate the $\vec{X}$ from a given $\vec{x}$. In our case, we have the 250 sample points from the cosine function so let's use that (the parameters and everything will be redefined and evaluated for clarity). Also, we are going to use the exponential representation for a change (nothing different as Euler's formula ensures that they are the same)

_Disclaimer: Even though $\{X_k\}$ have their values defined without the $\frac{1}{N}$, from here on, we'll just apply the division for practicality as if we're defining and using some new $X_k' \equiv \frac{X_k}{N}$ in place so we can write:_

$$x_n = \sum_{k=0}^{N-1} X_k' \cdot e^{\,i\frac{ 2 \pi}{N}kn}$$

dropping the $1/N$ factor outside of the summation symbol.

So we did the transformation and now we have a complex array $\vec X$ of the same size as $\vec x$ on our hands. So what? What did we accomplished?

Before delving into the meaning (or usefulness) of what we did, let's plot the real and imaginary parts of the array:

It seems like, major part of the components are shadowed by a couple of entries in comparison. Let's find out which entries they are:

As can be seen, the 12th and the 238th entries differ only by the sign of the imaginary component (i.e., they are each other's complex conjugates) and have the greatest magnitude. The 2nd highest values' (top entries' consecutive entries 13th and 237th) rapidly drop to 1/20th.

The listing looks very interesting as if there is a mirror in the center (indice = 125 : all the pair members have the same distance to this point, e.g., |12--125--238|, |13--125--237|, |11--125--239|, etc. ) that changes the sign of the imaginary part from one side to the other site of it.

As these are our $\{X_k\}$s, and the signal is the same one as the one we reproduced via "hands only" calculations above, let's think about what it means. We had -very succesfully- derived $(1.6 + 0.i)$ and $(1.6 - 0.i)$ for the $X_{12}$ and $X_{238}$ (we will later talk about it being 238 instead of -12, please be patient for now). However, here the dominant contributions seem to have the amplitudes of $1.579+i0.240$ and $1.579-i0.240$, respectively. Not exactly the same but close. As the remaining ones are much lower in amplitude, let's ignore them and try to reconstruct the signal using only the 12th and 238th entries:

As you can observe around the left and right boundaries, it is not an exact fit, albeit a good one. If we plot the errors:

Now, we'll deal with two issues:

  1. We used "real" parts of the reconstructed signal -- why did we feel the need? Wasn't it guaranteed that we'd get a real result?
    We did actually get a real result, alas since we were dealing with complex numbers, even if the imaginary part is zero (or something like zero) the number will be returned as complex. Just to make sure, checking the obtained signal's imaginary parts' magnitudes yields that they are practically zero indeed:
  1. Shouldn't we have gotten a perfect match?
    We would, if we were to include all the contributions:

And here's a follow-up question:

The Importance of Sampling

We have already talked about the importance of sampling in the beginning from the point of the quality of the observed data (discrete) wrt the real signal (continous). The more measurements we made, the more the representation was faithful to the original.

From the point of transformation, number of samples possess an additional importance: the number of $X_k$s is equal to the number of $x_k$s, but their domains are different - the more numbers we have, the more sharper the distribution of the coefficients will be (they are in fact, in a way, histograms) and also will include more contributions entering into the summation concerning the addition over $n$ for each $X_k$.

To see this, let's re-transform our signal of $x(t) = 3.2\cos(2\pi (3) t)$ but this time sample it with $N=1000$ points instead of $N=250$ points like we did last time:

(The above calculation may take longer than the previous calculations... but no worries, once we learn about the Fast Fourier Transfrom (FFT), they will be calculated in a jiff!)

And here are the top 10 contributions, sorted. Once again we have the 12th and this time 1000 - 12 = 988th coefficients.

We will return to the importance of sampling once more when we deal with a mixed signal (as opposed to a pure sinusoidal one we have been tackling with), to see that greater number of samples does also remove a faulty splitting in frequencies (but you didn't hear it from me - yet ;)

The Meaning of $\{X_k\}$ and the domain it lives

Let's reconstruct the signal and our coefficients from scratch, using $N=250$ samples as we did:

Real or Imaginary?

Here, the real parts' magnitudes are greater than the imaginary parts because they are related to cosine affinity whereas the imaginary coefficients are related to sines.

This is due to the fact that when we add up two corresponding frequencies' contributions, as our signal is real, the imaginary parts should cancel each other and given that the coefficients $X_k$ and $X_{N-k}$ are complex conjugates of each other in the form $a+ib$ and $a-ib$, we end up with:

$$\begin{gather*} \underbrace{(a+ib)}_{X_k}[\cos\theta + i\sin\theta] + \underbrace{(a-ib)}_{X_{N-k}}[\cos(-\theta) + i\sin(-\theta)]\\ =(a+ib)[\cos\theta + i\sin\theta] + (a-ib)[\cos\theta - i\sin\theta]\\ =(a\cos\theta+ia\sin\theta+ib\cos\theta-b\sin\theta)\\+ (a\cos\theta -ia\sin\theta-ib\cos\theta-b\sin\theta)\\ =\boxed{2a\cos(\theta) - 2b\sin(\theta)} \end{gather*}$$

For the pure cosine wave we have: $X_{12} = 1.579+i0.240$ and $X_{238} = 1.579-i0.240$ meaning $3.158\cos(\dots)-0.48\sin(\dots)$ Not the perfect result when compared to our actual signal of $3.2\cos(\dots) - 0.0\sin(\dots)$ but this is due to the relatively low number of samples. You can verify that with a higher number of samples, the coefficients converge to $(3.200, 0.000)$ (when we went from $N=250$ to $N=1000$ above, we had achieved $(3.198,0.120)$ already!).

And if we had a pure sine function, then the spikes would be maximized in the imaginary section. You can observe this by yourselves by running the following code (everything we've done so far, with only the "cos" replaced by "sin")

N = 250 # Number of samples
t = np.linspace(0,4,N)
T = 1/3 # Period
A = 3.2 # Amplitude

i2pi_N = 2j*np.pi/N

xs = A * np.sin(2 * np.pi / T * t)
Xs = np.zeros(N,complex)

for k in range(N):
    for n in range(N):
        Xs[k] += xs[n]*np.exp(-i2pi_N*k*n)

Xs /= N

plt.plot(np.real(Xs),"b-")
plt.plot(np.imag(Xs),"r-")
plt.legend(["Real Part","Imaginary Part"])
plt.show()

image-2.png

(the mysterious $(\theta)$ parameters will be revealed in the "One last remark" section once we get past the 'negative' frequencies)

'Take two'

Checking the formulation of the inverse transformation we see that $X_k$ plays the role of amplitude (& phase):

$$x_n = \frac{1}{N} \sum_{k=0}^{N-1} X_k\cdot e^{i 2 \pi k n / N}= \frac{1}{N} \sum_{k=0}^{N-1}X_k\cdot\left[\cos\left(\frac{2 \pi}{N}kn\right) + i \cdot \sin\left(\frac{2 \pi}{N}kn\right)\right]$$

As $|X_{12}|$ is much greater than all the remaining, let's assume that, this (along with the 238. component) is the only contributing factor and all the remaining coefficients are zero (/ignorable), and reconstruct $\vec x$ (as xx) using only these two components:

Our initial signal was, naturally, real. However from the formula, contribution from the imaginary part can be seen... but is it so? Let's check the constructed signal:

As can be seen, the imaginary parts are neglegible. With the highest value being on the order of $10^{-13}$:

We had some brief talk about this before but now we'll delve on it a bit more, so how did this happen? Remember the $X_k$ coefficients? They were consisting of complex conjugates, with $X_1 = X_{249}^*$, $X_2 = X_{248}^*$, ..., $X_{12} = X_{238}^*$, ... so let's see what happens when we multiply $X_{12}$ with $e^{-i{\frac{2\pi}{250}(12)n}}$ and $X_{238}$ with $e^{-i{\frac{2\pi}{250}(238)n}}$. Before we do that, let's check $X_{12}$ and $X_{238}$ explicitly:

so you see that, the complex conjugates runs to the rescue and gets us rid of the imaginary parts therefore no problems in ignoring them:

Let's compare this re-constructed $\vec {xx}$ signal with the original $\vec x$:

It's a good match! 8)

What about the horizontal axis’ meaning?

We still need to fix the horizontal axis' meaning. Remember that, after the Fourier transformation, the indices with the highest magnitude were the 12th and the 238th (out of 250 entries). We also observed that there was a symmetrical distribution: components with the same magnitudes appeared on the both sides of the center entry #125. This is not a coincidence.

Instead of adding the contributions from the 12th and 238th entries, let's construct separate signals from each, as if only that one exists. We'll label these as xx_12 and xx_238:

They are yielding exactly the same signal (half in amplitude of the original one).

What this means is: for a real signal, the first half of the coefficients carry all the information because, the other half will be complex conjugates of them and the imaginary parts will cancel each other out when reconstructing the signal (they ought to, as our original signal is real). Due to the periodicity, we can actually shift the right part to the left of the first part and associate them with "negative" frequencies (in reality, we are just changing our focus interval from $[0,N)$ to $[-N/2,N/2)$). As for the real signals, the magnitudes of their complex conjugate components are the same, we could just reconstruct the signal by using half of the interval's coefficients and then multiplying it by 2.

Moreover, looking at the plot on the left, we see that, the signal constructed from the $k=12$ (or $k=238$ for that matters) oscillates 12 times, not surprisingly, as $\tfrac{k}{N}$ in $\left(\frac{2\pi n k}{N}\right)$ acts directly as the frequency. But our original signal's frequency wasn't 12, it was set to 3 Hz. So there must be a normalization issue, which is due to the fact that when sampling, we sampled during 4 seconds, so there were 12 oscillations per 4 seconds resulting the frequency being 3 Hz (meaning "3 oscillations per second"). We didn't feed this information to Fourier transformation: we just supplied the 250 entries' displacements without specifying how much time it took. The calculated frequency is thus unnormalized -- we must divide the horizontal scale by 4 to have it correspond to the actual scale:

or if we shift the second half of the frequencies as negative frequencies (as our signal is of periodicity N, we can either take [0,𝑁)[0,N) or [−𝑁/2,𝑁/2)[−N/2,N/2) (or any other sequential 𝑁N element) intervals):

Zooming in:

There's actually a dedicated function in numpy (and/or scipy) to make this frequency adjustment calculation for us. It is the fftfreq() function and works by inputting the number of samples along with the spacing (/distance) between them. In our case we have 250 samples and there is a space of 4 / 250 between them (the spacing is equal to the inverse of the sampling rate: we have done 250 sampling in 4 seconds so our rate is 250/4 samples/second and its inverse is 4/250), so:

Therefore, the {X_k} coefficients are actually a measure of which frequency contributes how much and once we give the real meaning to the horizontal axis, we get the precious frequency spectrum of the signal.

As we have seen and derived and proved, for a real signal, no new information is contained in the 'negative' frequency range, therefore we can just multiply the positive ones by 2 and reconstruct the same signal. So instead of the above spectrum, we have the same information below:

How come we can shift the frequencies?

When we were talking about the horizontal axis' meaning and how to obtain the correct window with respect to the sampling rate, we easily shifted the frequencies on the right side to the negative region.

We've already saw that, when the signal is real, half of the frequencies' coefficients will be equal to the other half's coefficients' complex conjugate values. Furthermore, these corresponding pairs' ($k, N-k$) frequencies will be "inverse" to each other. (What does a 'negative frequency' mean by the way?! 8)

For the moment, let's show that waves with frequencies $\left(\frac{k}{N}\right)$ and $\left(\frac{k}{N} + mN\right)$ are equal to each other (with $\{k,N,m,n\}\in\mathbb{Z}$, $m$ arbitrary, $\{k,n\}=0\dots N-1$), so here we go:

$$\begin{align*}\cos\left[2\pi\left(\frac{k}{N} + mN\right) n\right] &= \cos\left(\frac{2\pi kn}{N} + 2\pi mNn\right)\\ &=\cos\left(\frac{2\pi kn}{N}\right)\cos\left(2\pi mNn\right)-\sin\left(\frac{2\pi kn}{N}\right)\sin\left(2\pi mNn\right) \end{align*} $$

where we have used the trigonometric identity: $\cos(\alpha+\beta) \equiv \cos\alpha \cos\beta - \sin\alpha \sin\beta$. Out critical term is the $\left(2\pi mNn\right)$ and more specifically $mNn$ in there. As all $m,N,n$ are integers, that parantheses is nothing but $2\pi$ multiplied by an integer, say $a$, and we all know how that goes for cosine and sine (given $a\in\mathbb{Z}$):

$$\cos(2\pi a) = 1,\quad \sin(2\pi a) = 0\quad\quad\left(a\in\mathbb{Z}\right)$$

Placing this equivalency in the above equation, we get: $$\begin{align*}\cos\left[2\pi\left(\frac{k}{N} + mN\right) n\right] &=\cos\left(\frac{2\pi kn}{N}\right)\underbrace{\cos\left(2\pi\overbrace{ mNn}^{\in\mathbb{Z}}\right)}_{1}-\sin\left(\frac{2\pi kn}{N}\right)\underbrace{\sin\left(2\pi \overbrace{ mNn}^{\in\mathbb{Z}}\right)}_{0}\\ &=\cos\left(\frac{2\pi kn}{N}\right)=\cos\left[2\pi\left(\frac{k}{N}\right)n\right] \end{align*}$$

Hence, we can add/subtract as many $N$ as we want to/from the frequency (thus shifting it) as long as we do it in integer multiplers of $N$.

When we look at it this way, maybe it isn't that impressive, however, from applicational point of view, we are boldly stating that even if the frequency is changed by a significant amount, no additional information will be gained regarding the data.

To see how this works, we will consider a very small system with $N = 5$ and $k=3$. The corresponding frequency will be: $\frac{k}{N}=\frac{3}{5}$, so let's plot it:

As we know the underlying frequency, let's also put it in:

Keep in mind that, we are making our calculations exactly at the integer values ($n=0,\dots,N-1$). The pale continous line (whose frequency is equal to $\tfrac{k}{N}$) is something we are not interested except at the integer parameter values.

Obviously, two waves with different frequencies will have very distinct graphs (duh!). As the frequency increases, the number of oscillations per unit length/time will also increase. So, let's plot the above wave of frequency $\tfrac{k}{N} = \tfrac{3}{5}$ along with another wave of frequency $\tfrac{k}{N}+mN = \tfrac{3}{5}+(2\times 5) = \tfrac{53}{5}$:

Nothing surprising here, as they are two waves as different as two waves with different frequencies can be! So how in the first place we could believe that there is no difference and we can take these as equivalent?!!!

The answer lies in the fact that we are only taking the integer parameters into account. Check the above graph with the points placed:

Zooming in:

meaning, as long as we are only interested in the integer parameters, all the $f = \tfrac{k}{N} + mN$ frequencies point to the same data, same information! (This fact is utilized beautifully in solid state physics, when we are dealing with phonons.)

One last remark..

To see the positive and "negative" frequencies in action, Let's start from the beginning. Let's take our sample case once again:

$$x(t) = 3.2 \cos(2\pi (3) t)$$

from which we see that: Amplitude (A): 3.2 Frequency (f): 3

Now, let's try to write it as an addition of two complex waves with their amplitudes being complex conjugates of each other and also their frequencies being negative of each other:

$$\begin{align*}y_1 &= \left(a+ib\right)\left[\cos\left(2\pi (f) t\right) + i \sin\left(2\pi (f) t\right)\right]\\ y_2 &= \left(a-ib\right)\left[\cos\left(2\pi(-f) t\right) + i \sin\left(2\pi (-f) t\right)\right] \end{align*}$$

Before we proceed further, please check for yourself to see that the amplitudes of the two complex waves are indeed complex conjugates of each other and also the frequencies are negative of each other as well.

Let's add these two waves. Expanding the parantheses:

$$\begin{gather*}y_1 + y_2 = \left(a+ib\right)\left[\cos\left(2\pi (f) t\right) + i \sin\left(2\pi (f) t\right)\right] + \left(a-ib\right)[\underbrace{\cos\left(2\pi(-f) t\right)}_{\cos\left(2\pi(f) t\right)} + i \,\underbrace{\sin\left(2\pi (-f) t\right)}_{-\sin\left(2\pi (f) t\right)}]\\ =\left[a\cos\left(2\pi f t\right)+ia\sin\left(2\pi f t\right)+ib\cos\left(2\pi f t\right)-b\sin\left(2\pi f t\right)\right] + \left[a\cos\left(2\pi f t\right)-ia\sin\left(2\pi f t\right)-ib\cos\left(2\pi f t\right)-b\sin\left(2\pi f t\right)\right]\\ =\boxed{2a\cos\left(2\pi f t\right) - 2b\sin\left(2\pi f t\right)} \end{gather*}$$

Comparing this with the above $x(t) = 3.2 \cos(2\pi (3) t)$, it is obvious that:

If our initial pure signal was that of a sine wave, say:

$$x(t) = 3.2 \sin(2\pi (3) t)$$

then, to comply, our positive and negative frequency pair would obviously be:

$$\begin{align*}y_1 &= (-1.6i)\left[\cos(2\pi(3)t) + i\sin(2\pi(3)t)\right]\\ y_2 &= (1.6i)\left[\cos(2\pi(-3)t) + i\sin(2\pi(-3)t)\right]\\&= (1.6i)\left[\cos(2\pi(3)t) - i\sin(2\pi(3)t)\right]\end{align*}$$

Fast Fourier Transform (FFT)

In the above "Manual Transform" section, we have computed the Fourier transform using the classic formulation. However, with the advance of computing power and development of sparse matrix toolkits, a faster way $\left(O\left(N^2\right)\Rightarrow O(N \log N)\right)$) for Fourier transform has been devised in the 1960s by James Cooley and John Tukey.

We won't be dealing with how the method works but will directly benefit from it. Using it is as simple as typing fft:

Case Studies

1. The brotherhood of cosines

Up until now, we had a pure sinusoidal (cosine) signal, let's pepper up the signal by combining 3 waves:

So, there we have it, exactly as given. However, we did it with a high number of samples ($N=48000$). Which brings us to the next lesson (again)...

The importance of the number of samples (in regard to frequency splittings)

The amplitudes are accurate when the number of samples are high. Let's quickly re-do the above processes this time with a lower N:

as can be seen there is a splitting of the high frequency. (677, but also 678)

Finding the amplitudes (only cosines case)

Now that we have the frequencies pinned down, it is more or less straightforward to fine tune the amplitudes via least squares fitting method (or any other optimization algorithm). We have 3 unknown amplitudes $\vec{A} = \{A^{(0)}, A^{(1)}, A^{(2)}\}$ for the three frequencies $\{f^{(0)}, f^{(1)}, f^{(2)}\}$ and $n$ equations such that:

$$\begin{align*}A^{(0)}\cos(2\pi f^{(0)}t_0) + A^{(1)}\cos(2\pi f^{(1)}t_0) + A^{(2)}\cos(2\pi f^{(2)}t_0) &= x_0\\ \vdots \\ A^{(0)}\cos(2\pi f^{(0)}t_i) + A^{(1)}\cos(2\pi f^{(1)}t_i) + A^{(2)}\cos(2\pi f^{(2)}t_i) &= x_i\\ \vdots \\ A^{(0)}\cos(2\pi f^{(0)}t_{n-1}) + A^{(1)}\cos(2\pi f^{(1)}t_{n-1}) + A^{(2)}\cos(2\pi f^{(2)}t_{n-1}) &= x_{n-1} \end{align*}$$

which can be written as a matrix multiplication:

$$\underbrace{\begin{pmatrix}\cos(2\pi f^{(0)}t_0) && \cos(2\pi f^{(1)}t_0) && \cos(2\pi f^{(2)}t_0)\\ \vdots&&\vdots&&\vdots\\ \cos(2\pi f^{(0)}t_{n-1}) && \cos(2\pi f^{(1)}t_{n-1}) && \cos(2\pi f^{(2)}t_{n-1})\end{pmatrix}_{(n\times3)}}_{\mathbb{C}} \underbrace{\begin{pmatrix}A^{(0)}\\A^{(1)}\\A^{(2)}\end{pmatrix}_{(3\times1)}}_{\vec{A}} = \underbrace{\begin{pmatrix}x_0\\\vdots \\x_i\\\vdots \\x_{n-1}\end{pmatrix}_{(n\times1)}}_{\vec{x}} $$$${\mathbb{C}}\cdot\vec{A}=\vec{x}$$

This is an overdetermined set of linear equations (3 unknowns and n>3 equations), so least squares fitting is a suitable choice of attack:

Mind though, this working out the amplitudes straightforward is not a general approach, it worked here because there were only cosines involved. But no worries, before the lecture ends, we'll have a general approach (involving least squares fitting) in our hands.

Moral of the story: Always eyeball the calculations, make your estimations before you see the results and don't leave it all to computer, because the computer always returns a result.

Secondary -scientific- moral of the story: The frequencies are the essential targets - once you find them, correct amplitudes can be derived easily, so, focus on the frequencies, and the rest will follow.

2. Noisy signal

When we filtered out the "weak" frequencies in the spectrum above (the xt_norm_low>0.25 line), we started something new: we based our reasoning that the signal was composed of a couple of "strong" signals and the rest was "noise". So, let's add some actual random noise to our signal:

Fourier transform reveals the active frequencies:

Job done!

3. A slightly more complex case (with sine added)

So far, we've dealt with a single, pure cosine function, now let's build a signal consisting of multiple functions with different amplitudes and frequencies:

So, this is our signal, in all its glory. Our "recording" duration is 2 seconds, with 4000 samples (we have forgotten its constituents: as of now we don't know the involved functions, their frequencies, their amplitudes, not even how many there are).

Zooming in, it looks like this:

So, let's apply FFT to the signal:

Pick the peaks from the positive range:

Instead of problem-specific filtering with amplitudes greater than 1000, a general solution of filtering above 50% of the normalized amplitudes can be used as well. In order to do so, we normalize the amplitudes and take the ones that are greater than 0.5 (i.e., > 50%) (We have already did it like this in the above cases, and this is the reason why we did it in that way: to generalize the code)

The frequencies are caught very well except for the smearing of the highest frequency. This kind of "splittings" may occur near the high frequency range with a low number of samples (as we've seen and investigated in the "The importance of the number of samples (in regard to frequency splittings)" section above).

We have the frequencies, have the amplitudes and also can reconstruct the signal using the inverse fourier transform (mind how we multiply the coefficients by 2 to also count for the negative frequency contributions)

as can be seen, the errors are significantly higher on the edges. One way to remedy -as we had demonstrated in the "The importance of the number of samples (in regard to frequency splittings)" section, is to increase the number of samples. To see this, just increase the $N$ parameter at the beginning of this section (say 10 fold) and re-calculate - here it is:

But what happens if we do not have the opportunity to increase the number of samples (e.g., when we have already made our observations and it would be too impractical to re-do the measurements)? Then, we "use our heads" 8)

Optimizing the amplitudes once we have obtained the frequencies (the real deal)

Let's go once again to our relatively low number of samples and take it from there: Remember that we had filtered out the frequencies with normalized amplitudes less than 0.5? Let's decrease that threshold to 0.1 and then try to reconstruct the original signal:

Really much better indeed! Although by including all the "side" frequencies, we are actually losing the simplicity of the real signal. Even though the mathematical result is better, the physical point of view got much complicated. (This is the part where we'll start using our heads! 8)

Let's remember the frequencies and their corresponding amplitudes when we went for 50% amplitude filtering:

The amplitudes are not that well spotted but now that we have the frequencies, they can be very easily be fitted (using a method like least squares fit - in a minute).

Speaking of frequencies, we see that the obtained frequencies of 900 and 900.5 are too close to be a coincidence, therefore we'll keep the 900 and leave 900.5. This leaves us with 4 frequencies (and even though we "don't know" it (as we had made ourselves forgotten once we built the original signal from them), they are the exact frequencies we had used). We have seen that, with the calculated amplitudes they yielded not-so-great results, therefore, we'll play with the amplitudes using the least-squares fitting method:

At the moment, our summed up function looks like this:

and how different (/erroneous) is it from the original signal? Let's find out, but for practical purposes, let's define a function that takes these amplitudes as parameters (in the form of an array) and returns the difference:

or, alternatively (& equivalently), using the exponential form:

and here it is in action:

A0 is the array (/vector) that holds the amplitudes we had from the FFT, those will play a role as our starting point. We have included the real and imaginary parts as separate entries as most of the optimization functions are not supporting direct processing of complex parameters.

Our aim is to adjust the amplitudes vector to yield the minimum return (error) from our difference function. For this purpose we will use the least squares fitting method:

Getting rid of zeros (<1E-10):

Compare these with the input parameters:

Negative coefficient for the sine???

At this moment, you might be wondering why did we implemented the sin part to the function as
sig_con -= amp_sin*np.sin(2*np.pi*freqs_picked[i]*t) and also -directly as a result of that- got the amplitude as -1.13, instead of our input of 1.13. The reason lies in the expansion of the amplitudes in the complex form - the details are given and derived in the "One Last Remark" section above, but the summary is the equivalency of:

$$\begin{align*}y_1 &= \left(a+ib\right)\left[\cos\left(2\pi (f) t\right) + i \sin\left(2\pi (f) t\right)\right]\\ y_2 &= \left(a-ib\right)\left[\cos\left(2\pi(-f) t\right) + i \sin\left(2\pi (-f) t\right)\right] \end{align*}$$$$\begin{gather*}y_1 + y_2 = \left(a+ib\right)\left[\cos\left(2\pi (f) t\right) + i \sin\left(2\pi (f) t\right)\right] + \left(a-ib\right)[\underbrace{\cos\left(2\pi(-f) t\right)}_{\cos\left(2\pi(f) t\right)} + i \,\underbrace{\sin\left(2\pi (-f) t\right)}_{-\sin\left(2\pi (f) t\right)}]\\ =\left[a\cos\left(2\pi f t\right)+ia\sin\left(2\pi f t\right)+ib\cos\left(2\pi f t\right)-b\sin\left(2\pi f t\right)\right] + \left[a\cos\left(2\pi f t\right)-ia\sin\left(2\pi f t\right)-ib\cos\left(2\pi f t\right)-b\sin\left(2\pi f t\right)\right]\\ =\boxed{2a\cos\left(2\pi f t\right) - 2b\sin\left(2\pi f t\right)} \end{gather*}$$

so, the amplitudes of the sine components come as subtraction! ;) (Mind that, when we have a signal made of pure cosines and sines, the magnitude of the amplitudes are the absolute values of the coefficients, hence the positive amplitudes: a complex coefficient like $0-1.13i$ has a magnitude of $1.13$). To see this more clearly, imagine that we only had a sine wave as the signal, say $1.3\sin(2\pi(576)t)$: Fourier transform would pick positive and negative frequencies in the form of:

$$(a + ib)[\cos(2\pi(576)t) + i\sin(2\pi(576)t)] \\+ (a - ib)[\cos(2\pi(-576)t) + i \sin(2\pi(-576)t)] $$

and we would equate this to our input signal of $1.3\sin(2\pi(576)t)$. Expanding the complex coefficient over the $\cos(\dots) + i \sin(\dots)$ and also using $\cos(-\theta)=\cos(\theta),\;\sin(-\theta)=-\sin(\theta)$ the addition of the positive and negative frequencies contributions would be:

$$2a\cos(2\pi(576)t) - 2b\sin(2\pi(576)t)$$

and equating it to our original wave, we'd have:

$$2a\cos(2\pi(576)t) - 2b\sin(2\pi(576)t) = 1.3\sin(2\pi(576)t)$$

From here, we see that $a=0$ while $-2b = 1.3$ hence the coefficient $2b$ is found out to be: $2b = -1.3$.

Let's rebuilt the signal from the obtained ones:

and there we have it! <dance! dance! dance!>