Fast Fourier Transform (FFT) Algorithm Implementation In Python
Fast Fourier Transform (FFT) are used in digital signal processing and training models used in Convolutional Neural Networks (CNN). We demonstrate how to apply the algorithm using Python.
FFT has many applications in component analysis.
One application of FFT is in digital signal processing. It can be used as a tool for decomposing signals from the time domain to the frequency domain. This is a good way of separating signals like when quantizing analog waves to discrete signals.
Another application of FFT is in image processing. In convolutional layer overlays a kernel can scan a section of an image. This takes all values from that location to perform operations. The kernel then shifts to another section of the image and repeats the process until it scans the entire image.
Discrete Fourier Transforms (DFT)
FFT is also the computation of Discrete Fourier Transforms (DFT) into an algorithmic format. It is a way to compute Fourier Transforms much faster than using the conventional method. Let us take a look at DFT to get some idea of how FFT makes computing results faster.
DFT transforms a sequence of N numbers …
… to a sequence of another set of complex numbers:
Without going over the entire theorem, DFT is basically taking any quantity or signal that varies over time and decomposing it into its components (e.g. frequency). This can be the pixels in an image or the pressure of sound in radio waves (used in DSP chips).
The problem with DFT is that it takes more calculations to arrive at an answer. This is because you have to do N(multiplications) with N(additions), which in Big O Notation is O(N²) in terms of time complexity.
Fast Fourier Transforms (FFT)
With FFT we save some steps by reducing the number of calculations in DFT. As a result we can reduce a domain of size N from O(N²) to O(Nlog2N).
DFT using FFT can be written using the following formula:
We have the discrete signal x(n) multiplied with e (raised to a function specified) , with N representing the size of the domain for the results of the sum of a value n.
As the name implies, it is a fast or a faster way to compute Fourier Transforms. This becomes more noticeable when the size of the domain becomes larger. Take for example this table (taken from Towards Data Science article by Cory Maklin):
Let us say it takes 1 nanosecond per operation performed in a system. With FFT, it takes approximately 30 seconds to complete operations of size N = 10⁹. If FFT was not used, it would take 31.7 years (10¹⁸) using the conventional algorithm formula.
N = 10⁹ = 1 000 000 000 operationsWithout FFTN² = 10¹⁸ = 1 000 000 000 000 000 000 / 1 000 000 000 = 31.7 yearsNlog2N = 30 x 1 000 000 000 / 1 000 000 000 = 30 seconds
Implementing The Algorithm
Let us start with a simple example. From a Python environment at the prompt (you can also write a Python or py file), import the numpy library.
import numpy as npWe will now specify a function and use matrix operations to arrive at the result of the computation.
def dft(x):
x = np.asarray(x, dtype=float)
N = x.shape[0]
n = np.arange(N)
k = n.reshape((N, 1))
M = np.exp(-2j * np.pi * k * n / N)
return np.dot(M, x)Now let us generate a random array of numbers to use with the calculations.
x = np.random.random(1024)Check the value of x:
xarray([0.2100223 , 0.76314102, 0.45883551, ..., 0.75090954, 0.01397708, 0.66781247])
The next example makes use of a notebook using a recursive approach.
import matplotlib.pyplot as plt
import numpy as np
plt.style.use('seaborn-poster')
%matplotlib inlinedef FFT(x):
"""
A recursive implementation of
the 1D Cooley-Tukey FFT, the
input should have a length of
power of 2.
"""
N = len(x)
if N == 1:
return x
else:
X_even = FFT(x[::2])
X_odd = FFT(x[1::2])
factor = \
np.exp(-2j*np.pi*np.arange(N)/ N)
X = np.concatenate(\
[X_even+factor[:int(N/2)]*X_odd,
X_even+factor[int(N/2):]*X_odd])
return X# sampling rate
sr = 128
# sampling interval
ts = 1.0/sr
t = np.arange(0,1,ts)
freq = 1.
x = 3*np.sin(2*np.pi*freq*t)
freq = 4
x += np.sin(2*np.pi*freq*t)
freq = 7
x += 0.5* np.sin(2*np.pi*freq*t)
plt.figure(figsize = (8, 6))
plt.plot(t, x, 'r')
plt.ylabel('Amplitude')
plt.show()
This produces the graph plot:
Next we can calculate the Fourier Transform of the signal from the graph.
X=FFT(x)
# calculate the frequency
N = len(X)
n = np.arange(N)
T = N/sr
freq = n/T
plt.figure(figsize = (12, 6))
plt.subplot(121)
plt.stem(freq, abs(X), 'b', \
markerfmt=" ", basefmt="-b")
plt.xlabel('Freq (Hz)')
plt.ylabel('FFT Amplitude |X(freq)|')
# Get the one-sided specturm
n_oneside = N//2
# get the one side frequency
f_oneside = freq[:n_oneside]
# normalize the amplitude
X_oneside =X[:n_oneside]/n_oneside
plt.subplot(122)
plt.stem(f_oneside, abs(X_oneside), 'b', \
markerfmt=" ", basefmt="-b")
plt.xlabel('Freq (Hz)')
plt.ylabel('Normalized FFT Amplitude |X(freq)|')
plt.tight_layout()
plt.show()
FFT Implementation In Python
Here are more examples using Python libraries, with numpy and scipy. We used numpy earlier, but here are more advanced examples.
Let us first generate the signal.
import matplotlib.pyplot as plt
import numpy as np
plt.style.use('seaborn-poster')
%matplotlib inline# sampling rate
sr = 2000
# sampling interval
ts = 1.0/sr
t = np.arange(0,1,ts)
freq = 1.
x = 3*np.sin(2*np.pi*freq*t)
freq = 4
x += np.sin(2*np.pi*freq*t)
freq = 7
x += 0.5* np.sin(2*np.pi*freq*t)
plt.figure(figsize = (8, 6))
plt.plot(t, x, 'r')
plt.ylabel('Amplitude')
plt.show()
The signal is identical to the previous recursive example.
We will now use the fft and ifft functions from numpy to calculate the FFT amplitude spectrum and inverse FFT to obtain the original signal. We then plot both results and time the fft function using a 2000 length signal.
from numpy.fft import fft, ifft
X = fft(x)
N = len(X)
n = np.arange(N)
T = N/sr
freq = n/T
plt.figure(figsize = (12, 6))
plt.subplot(121)
plt.stem(freq, np.abs(X), 'b', \
markerfmt=" ", basefmt="-b")
plt.xlabel('Freq (Hz)')
plt.ylabel('FFT Amplitude |X(freq)|')
plt.xlim(0, 10)
plt.subplot(122)
plt.plot(t, ifft(X), 'r')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.tight_layout()
plt.show()
%timeit fft(x)We get the result:
36.2 µs ± 775 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)Next, we will use scipy instead of numpy.
from scipy.fftpack import fft, ifft
X = fft(x)
plt.figure(figsize = (12, 6))
plt.subplot(121)
plt.stem(freq, np.abs(X), 'b', \
markerfmt=" ", basefmt="-b")
plt.xlabel('Freq (Hz)')
plt.ylabel('FFT Amplitude |X(freq)|')
plt.xlim(0, 10)
plt.subplot(122)
plt.plot(t, ifft(X), 'r')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.tight_layout()
plt.show()
%timeit fft(x)We get the result:
14.8 µs ± 471 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)Synopsis
Using the FFT algorithm is a faster way to get DFT calculations. The built-in Python functions for FFT are quite fast and easy to use, notably the scipy library. It significantly lessens the amount of time, which can also save costs. These functions are available in Python, so developers can build ready-to-use systems that can utilize these library modules.
