We've all seen those weird looking mathematics equations that pop up when we hear about **Fourier transforms**. In this blog post, we'll programatically try and develop an intuitive understanding into the whole process. While I'll be using the scientific Python stack in this blog post, code in Matlab, R should not be that different.

First, let us assume that we are doing some signal acquisition and we can sample at 100 Hz frequency (100 times per second). We collect data for 10 seconds. So, we have a total of 1000 samples.

```
Samples collection duration (T) = 10s
Sampling frequency (Fs) =100Hz
Number of samples (N) = Ts*Fs=1000
```

It must be noted that Fs is the rate at which we are able to sample. Let us now create a signal composed of different frequencies.

In [1]:

```
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('fivethirtyeight')
import warnings
warnings.filterwarnings('ignore')
%matplotlib inline
```

In [2]:

```
def create_sine(t, fs, frequency_of_sine):
return np.sin(2*np.pi*frequency_of_sine*np.arange(0, t, 1./fs))
def create_cosine(t, fs, frequency_of_cosine):
return np.cos(2*np.pi*frequency_of_cosine*np.arange(0, t, 1./fs))
```

In [3]:

```
plt.figure(figsize=(8,4))
plt.plot(create_sine(10, 100, 1))
plt.xlabel("Samples")
plt.ylabel("Amplitude");
```

In [4]:

```
t = 10
fs = 100
N = t*fs
num_components = 4
components = np.zeros((num_components, N))
components[0] = np.ones(N)
components[1] = create_sine(t, fs, 10)
components[2] = create_sine(t, fs, 2)
components[3] = create_sine(t, fs, 0.5)
```

Plotting these four components individually.

In [5]:

```
fig, ax = plt.subplots(nrows=num_components, sharex=True, figsize=(12,6))
for i in range(num_components):
ax[i].plot(components[i])
ax[i].set_ylim((-1.1, 1.1))
ax[i].set_title('Component {}'.format(i))
ax[i].set_ylabel("Amplitude")
ax[num_components-1].set_xlabel("Samples")
plt.tight_layout()
```

Let our signal be

$$x = -0.5\times x_0 + 0.1\times x_1 + 0.2\times x_2 -0.6\times x_3 $$In [6]:

```
x = -0.5*components[0]+0.1*components[1]+0.2*components[2]-0.6*components[3]
```

Our dummy signal looks something like the following

In [7]:

```
plt.plot(x)
plt.xlabel("Samples")
plt.ylabel("Amplitude");
```

For finding the various frequency components in the signal, we'll be using the Discrete Fourier Transform (DFT). The key step in DFT is to find the `correlation`

between cosine waves of different frequencies with the signal that we intend to process. A high amplitude of this `correlation`

indicates the presence of this frequency in our signal.

It must be noted that the definition of `correlation`

is different from the definition we encounter in statistics. Here, `correlation`

between two signals simply means the dot product between the two. Please note that the dot product is the sum of the element wise product between the two signals.

Let us see the `correlation`

between a few pair of signals.

In [8]:

```
signal_1 = create_sine(10, 100,0.1)
signal_2 = create_cosine(10, 100,0.1 )
plt.plot(signal_1, label="Sine")
plt.plot(signal_2, label="Cosine")
plt.title("Correlation={}".format(np.sum(signal_1*signal_2)))
plt.legend();
```

In [9]:

```
signal_1 = create_sine(10, 100,0.1)
signal_2 = create_sine(10, 100,0.1)
plt.plot(signal_1, label="Sine")
plt.plot(signal_2, label="Sine 2",linestyle='--')
plt.title("Correlation={}".format(np.sum(signal_1*signal_2)))
plt.legend();
```

`correlation`

between our signal and a cosine wave of different frequencies. We also create small helper functions to create sine and cosine waves containing `k`

periods in `N`

sample points.

In [10]:

```
def create_cosine_k_N(k, N):
return np.cos((2*np.pi*k/N)*np.arange(N))
def create_sine_k_N(k, N):
return np.sin((2*np.pi*k/N)*np.arange(N))
```

`correlation`

between this cosine and our signal `x`

.

In [11]:

```
cos_5 = create_cosine_k_N(5,N)
plt.plot(x, label="x")
plt.plot(cos_5, label="cosine (0.5Hz)")
plt.plot(cos_5*x, label="cosine (0.5 Hz)*x")
plt.title("Correlation={}".format(np.sum(cos_5*x)))
plt.legend();
```

`correlation`

with the sine of the same frequency?

In [12]:

```
sin_5 = create_sine_k_N(5,N)
plt.plot(x, label="x")
plt.plot(sin_5, label="sine (0.5Hz)")
plt.plot(sin_5*x, label="sine (0.5 Hz)*x")
plt.title("Correlation={}".format(np.sum(sin_5*x)))
plt.legend();
```

Great! The correlation with a 0.5 Hz sinusoidal has a high amplitude. We're on to something now! Yes, we want to detect the presence of a particular frequency and not really worry about phase for the moment.

So, for each frequency that we find our correlation with the original signal, we can have two components- the correlation with the cosine (called the real component) and the correlation with the sine of that frequency (called the imaginary component). This is exactly what those weird looking DFT equations mean! A high absolute value of either of these components suggests the presence of that particular frequency in our signal.

So, for every frequency in 0 to N-1 Hz, we repeat this procedure and obtain the DFT coefficients. Finally, the mathematics looks something like the following:

$ \mathrm{DFT(K)} = \sum{x_n\times \cos{\frac{2\pi K}{N}}} - \iota \sum{x_n\times \sin{\frac{2\pi K}{N}}} $

Now, if you are guessing that there is some relation with Euler's formula, you are right!

In [13]:

```
plt.figure(figsize=(12,4))
sin_20 = create_sine_k_N(20,N)
plt.plot(x, label="x")
plt.plot(sin_20, label="sine (2 Hz)")
plt.plot(sin_20*x, label="sine (2 Hz)*x")
plt.title("Correlation={}".format(np.sum(sin_20*x)))
plt.legend();
```

How about the DC component?

In [14]:

```
plt.figure(figsize=(12,4))
cos_0 = create_cosine_k_N(0,N)
plt.plot(x, label="x")
plt.plot(cos_0, label="cosine 0 Hz)")
plt.plot(cos_0*x, label="cosine (0 Hz)*x")
plt.title("Correlation={}".format(np.sum(cos_0*x)))
plt.legend();
```

In [15]:

```
fft_x = np.fft.fft(x)
```

`fft`

routine returns an array of length 1000, which is equal to the number of samples. This refers to the 1000 DFT coefficients, from K=0 to 999.

In [16]:

```
len(fft_x)
```

Out[16]:

In [17]:

```
fft_x[0]
```

Out[17]:

The above was the DFT coefficient for K=0. It has two components (real and imaginary) arising from the cosine and the sine waves respectively. However, we are concerned with the amplitude of the signal. This can be done by considering the absolute value of these coefficients.

We will now plot these DFT coefficients.

In [18]:

```
plt.plot(abs(fft_x))
plt.xlim((-100, 1100))
plt.ylim((-5, 520))
plt.xlabel("K")
plt.ylabel("|DFT(K)|");
```

In [19]:

```
plt.plot(abs(fft_x)[:500])
plt.xlim((-5, 500))
plt.ylim((-5, 520))
plt.xlabel("K")
plt.ylabel("|DFT(K)|");
```

In [20]:

```
plt.plot(abs(fft_x)[:500])
plt.xlim((-5, 120))
plt.ylim((-5, 520))
plt.xlabel("K")
plt.ylabel("|DFT(K)|");
```

In [21]:

```
plt.plot(np.arange(0, 500)/10.,abs(fft_x)[:500])
plt.xlim((-1, 12))
plt.ylim((-5, 520))
plt.xlabel("Freq(Hz)")
plt.ylabel("|DFT(K)|");
```

In [27]:

```
Et = np.sum(x*x)
print Et
Ef = np.sum(abs(fft_x)*abs(fft_x))/N
print Ef
```

Energy of DC component

In [30]:

```
Et_DC = np.sum(-0.5*components[0]*-0.5*components[0])
ft_DC = np.fft.fft(-0.5*components[0])
Ef_DC = np.sum(abs(ft_DC)*abs(ft_DC))/N
print(Et_DC, Ef_DC)
```

Energy of 10 Hz component

In [31]:

```
Et_10hz = np.sum(0.1*components[1]*0.1*components[1])
ft_10hz = np.fft.fft(0.1*components[1])
Ef_10hz = np.sum(abs(ft_10hz)*abs(ft_10hz))/N
print(Et_10hz, Ef_10hz)
```

In [32]:

```
sqr_wave = np.hstack([np.ones(250), np.zeros(250), np.ones(250), np.zeros(250)])
```

In [33]:

```
plt.plot(sqr_wave)
plt.ylim((-0.2,1.2));
```

In [34]:

```
dft_sqr = np.fft.fft(sqr_wave)
```

In [35]:

```
plt.plot(np.arange(0, 500)/10.,abs(dft_sqr)[:500])
plt.xlim((-1, 12))
plt.ylim((-5, 520))
plt.xlabel("Freq(Hz)")
plt.ylabel("|DFT(K)|");
```

There you go! Now, does FT sound all that difficult?

I have used a number of wonderful references for developing my understanding. I'll try and list a few here: