Fourier transforms

Functions related to Fourier transforms can be called by importing the fft sub-module first.

numpy: https://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.ifft.html

fft

Since ulab’s ndarray does not support complex numbers, the invocation of the Fourier transform differs from that in numpy. In numpy, you can simply pass an array or iterable to the function, and it will be treated as a complex array:

# code to be run in CPython

fft.fft([1, 2, 3, 4, 1, 2, 3, 4])
array([20.+0.j,  0.+0.j, -4.+4.j,  0.+0.j, -4.+0.j,  0.+0.j, -4.-4.j,
        0.+0.j])

WARNING: The array that is returned is also complex, i.e., the real and imaginary components are cast together. In ulab, the real and imaginary parts are treated separately: you have to pass two ndarrays to the function, although, the second argument is optional, in which case the imaginary part is assumed to be zero.

WARNING: The function, as opposed to numpy, returns a 2-tuple, whose elements are two ndarrays, holding the real and imaginary parts of the transform separately.

# code to be run in micropython

import ulab as np
from ulab import vector
from ulab import fft

x = np.linspace(0, 10, num=1024)
y = vector.sin(x)
z = np.zeros(len(x))

a, b = fft.fft(x)
print('real part:\t', a)
print('\nimaginary part:\t', b)

c, d = fft.fft(x, z)
print('\nreal part:\t', c)
print('\nimaginary part:\t', d)
real part:   array([5119.996, -5.004663, -5.004798, ..., -5.005482, -5.005643, -5.006577], dtype=float)

imaginary part:      array([0.0, 1631.333, 815.659, ..., -543.764, -815.6588, -1631.333], dtype=float)

real part:   array([5119.996, -5.004663, -5.004798, ..., -5.005482, -5.005643, -5.006577], dtype=float)

imaginary part:      array([0.0, 1631.333, 815.659, ..., -543.764, -815.6588, -1631.333], dtype=float)

ifft

The above-mentioned rules apply to the inverse Fourier transform. The inverse is also normalised by N, the number of elements, as is customary in numpy. With the normalisation, we can ascertain that the inverse of the transform is equal to the original array.

# code to be run in micropython

import ulab as np
from ulab import vector
from ulab import fft

x = np.linspace(0, 10, num=1024)
y = vector.sin(x)

a, b = fft.fft(y)

print('original vector:\t', y)

y, z = fft.ifft(a, b)
# the real part should be equal to y
print('\nreal part of inverse:\t', y)
# the imaginary part should be equal to zero
print('\nimaginary part of inverse:\t', z)
original vector:     array([0.0, 0.009775016, 0.0195491, ..., -0.5275068, -0.5357859, -0.5440139], dtype=float)

real part of inverse:        array([-2.980232e-08, 0.0097754, 0.0195494, ..., -0.5275064, -0.5357857, -0.5440133], dtype=float)

imaginary part of inverse:   array([-2.980232e-08, -1.451171e-07, 3.693752e-08, ..., 6.44871e-08, 9.34986e-08, 2.18336e-07], dtype=float)

Note that unlike in numpy, the length of the array on which the Fourier transform is carried out must be a power of 2. If this is not the case, the function raises a ValueError exception.

spectrogram

In addition to the Fourier transform and its inverse, ulab also sports a function called spectrogram, which returns the absolute value of the Fourier transform. This could be used to find the dominant spectral component in a time series. The arguments are treated in the same way as in fft, and ifft.

# code to be run in micropython

import ulab as np
from ulab import vector
from ulab import fft

x = np.linspace(0, 10, num=1024)
y = vector.sin(x)

a = fft.spectrogram(y)

print('original vector:\t', y)
print('\nspectrum:\t', a)
original vector:     array([0.0, 0.009775015390171337, 0.01954909674625918, ..., -0.5275140569487312, -0.5357931822978732, -0.5440211108893639], dtype=float)

spectrum:    array([187.8635087634579, 315.3112063607119, 347.8814873399374, ..., 84.45888934298905, 347.8814873399374, 315.3112063607118], dtype=float)

As such, spectrogram is really just a shorthand for np.sqrt(a*a + b*b):

# code to be run in micropython

import ulab as np
from ulab import fft
from ulab import vector

x = np.linspace(0, 10, num=1024)
y = vector.sin(x)

a, b = fft.fft(y)

print('\nspectrum calculated the hard way:\t', vector.sqrt(a*a + b*b))

a = fft.spectrogram(y)

print('\nspectrum calculated the lazy way:\t', a)
spectrum calculated the hard way:    array([187.8641, 315.3125, 347.8804, ..., 84.4587, 347.8803, 315.3124], dtype=float)

spectrum calculated the lazy way:    array([187.8641, 315.3125, 347.8804, ..., 84.4587, 347.8803, 315.3124], dtype=float)

Computation and storage costs

RAM

The FFT routine of ulab calculates the transform in place. This means that beyond reserving space for the two ndarrays that will be returned (the computation uses these two as intermediate storage space), only a handful of temporary variables, all floats or 32-bit integers, are required.

Speed of FFTs

A comment on the speed: a 1024-point transform implemented in python would cost around 90 ms, and 13 ms in assembly, if the code runs on the pyboard, v.1.1. You can gain a factor of four by moving to the D series https://github.com/peterhinch/micropython-fourier/blob/master/README.md#8-performance.

# code to be run in micropython

import ulab as np
from ulab import vector
from ulab import fft

x = np.linspace(0, 10, num=1024)
y = vector.sin(x)

@timeit
def np_fft(y):
    return fft.fft(y)

a, b = np_fft(y)
execution time:  1985  us

The C implementation runs in less than 2 ms on the pyboard (we have just measured that), and has been reported to run in under 0.8 ms on the D series board. That is an improvement of at least a factor of four.