How to compute Discrete Fourier Transform (DFT) using SciPy

[box type=”note” align=”” class=”” width=””]This article is an excerpt from a book co-authored by L. Felipe Martins, Ruben Oliva Ramos and V Kishore Ayyadevara titled SciPy Recipes. This book provides numerous recipes to tackle day-to-day challenges associated with scientific computing and data manipulation using SciPy stack.[/box]

Today, we will compute Discrete Fourier Transform (DFT) and inverse DFT using SciPy stack. In this article, we will focus majorly on the syntax and the application of DFT in SciPy assuming you are well versed with the mathematics of this concept.

Discrete Fourier Transforms  

A discrete Fourier transform transforms any signal from its time/space domain into a related signal in frequency domain. This allows us to not only analyze the different frequencies of the data, but also enables faster filtering operations, when used properly. It is possible to turn a signal in a frequency domain back to its time/spatial domain, thanks to inverse Fourier transform (IFT).

How to do it…

To follow with the example, we need to continue with the following steps:

1. The basic routines in the scipy.fftpack module compute the DFT and its inverse, for discrete signals in any dimension—fft, ifft (one dimension), fft2, ifft2 (two dimensions), and fftn, ifftn (any number of dimensions).
2. Verify all these routines assume that the data is complex valued. If we know beforehand that a particular dataset is actually real-valued, and should offer realvalued frequencies, we use rfft and irfft instead, for a faster algorithm.
3. In order to complete with this, these routines are designed so that composition with their inverses always yields the identity.
4. The syntax is the same in all cases, as follows:

fft(x[, n, axis, overwrite_x])

The first parameter, x, is always the signal in any array-like form. Note that fft performs one-dimensional transforms. This means that if x happens to be two-dimensional, for example, fft will output another two-dimensional array, where each row is the transform of each row of the original. We can use columns instead, with the optional axis parameter. The rest of the parameters are also optional; n indicates the length of the transform and overwrite_x gets rid of the original data to save memory and resources. We usually play with the n integer when we need to pad the signal with zeros or truncate it. For a higher dimension, n is substituted by shape (a tuple) and axis by axes (another tuple). To better understand the output, it is often useful to shift the zero frequencies to the center of the output arrays with ifftshift. The inverse of this operation, ifftshift, is also included in the module.

How it works…

The following code shows some of these routines in action when applied to a checkerboard:

import numpy

from scipy.fftpack import fft,fft2, fftshift

import matplotlib.pyplot as plt

B=numpy.ones((4,4)); W=numpy.zeros((4,4))

signal = numpy.bmat("B,W;W,B")

onedimfft = fft(signal,n=16)

twodimfft = fft2(signal,shape=(16,16))

plt.figure()

plt.gray()

plt.subplot(121,aspect='equal')

plt.pcolormesh(onedimfft.real)

plt.colorbar(orientation='horizontal')

plt.subplot(122,aspect='equal')

plt.pcolormesh(fftshift(twodimfft.real))

plt.colorbar(orientation='horizontal')

plt.show()

Note how the first four rows of the one-dimensional transform are equal (and so are the last four), while the two-dimensional transform (once shifted) presents a peak at the origin and nice symmetries in the frequency domain.

In the following screenshot, which has been obtained from the previous code, the image on the left is the fft and the one on the right is the fft2 of a 2 x 2 checkerboard signal:

Computing the discrete Fourier transform (DFT) of a data series using the FFT Algorithm

In this section, we will see how to compute the discrete Fourier transform and some of its Applications.

How to do it…

In the following table, we will see the parameters to create a data series using the FFT algorithm:

How it works…

This code represents computing an FFT discrete Fourier in the main part:

np.fft.fft(np.exp(2j * np.pi * np.arange(8) / 8))

array([ -3.44505240e-16 +1.14383329e-17j,

8.00000000e+00 -5.71092652e-15j,

2.33482938e-16 +1.22460635e-16j,

1.64863782e-15 +1.77635684e-15j,

9.95839695e-17 +2.33482938e-16j,

0.00000000e+00 +1.66837030e-15j,

1.14383329e-17 +1.22460635e-16j,

-1.64863782e-15 +1.77635684e-15j])

In this example, real input has an FFT that is Hermitian, that is, symmetric in the real part and anti-symmetric in the imaginary part, as described in the numpy.fft documentation.

import matplotlib.pyplot as plt

t = np.arange(256)

sp = np.fft.fft(np.sin(t))

freq = np.fft.fftfreq(t.shape[-1])

plt.plot(freq, sp.real, freq, sp.imag)

[<matplotlib.lines.Line2D object at 0x...>, <matplotlib.lines.Line2D object

at 0x...>]

plt.show()

The following screenshot shows how we represent the results:

Computing the inverse DFT of a data series

In this section, we will learn how to compute the inverse DFT of a data series.

How to do it…

In this section we will see how to compute the inverse Fourier transform.

The returned complex array contains y(0), y(1),…, y(n-1) where:

How it works…

In this part, we represent the calculous of the DFT:

np.fft.ifft([0, 4, 0, 0])

array([ 1.+0.j, 0.+1.j, -1.+0.j, 0.-1.j])

Create and plot a band-limited signal with random phases:

import matplotlib.pyplot as plt

t = np.arange(400)

n = np.zeros((400,), dtype=complex)

n[40:60] = np.exp(1j*np.random.uniform(0, 2*np.pi, (20,)))

s = np.fft.ifft(n)

plt.plot(t, s.real, 'b-', t, s.imag, 'r--')

plt.legend(('real', 'imaginary'))

plt.show()

Then we represent it, as shown in the following screenshot:  

We successfully explored how to transform signals from time or space domain into frequency domain and vice-versa, allowing you to analyze frequencies in detail.

If you found this tutorial useful, do check out the book SciPy Recipes to get hands-on recipes to perform various data science tasks with ease.