Fast Fourier Transform

Mathematical Definitions

Fast Fourier Transforms are efficient algorithms for calculating the discrete Fourier transform (DFT)

\[x_j = \sum_{k=0}^{n-1} z_k \exp(-2\pi i j k / n)\]

The DFT usually arises as an approximation to the continuous Fourier transform when functions are sampled at discrete intervals in space or time. The naive evaluation of the discrete Fourier transform is a matrix-vector multiplication \(W\vec{z}\). A general matrix-vector multiplication takes O(n²) operations for n data-points. Fast Fourier transform algorithms use a divide-and-conquer strategy to factorize the matrix W into smaller sub-matrices, corresponding to the integer factors of the length n. If n can be factorized into a product of integers f₁ f₂ ... f_m then the DFT can be computed in O(\(n \sum f_i\)) operations. For a radix-2 FFT this gives an operation count of O(n log₂ n).

All the FFT functions offer two types of transform: forwards and inverse, based on the same mathematical definitions. The definition of the forward Fourier transform, fft(z), is,

\[x_j = \sum_{k=0}^{n-1} z_k \exp(-2\pi i j k / n)\]

and the definition of the inverse Fourier transform is,

\[z_j = {1 \over n} \sum_{k=0}^{n-1} x_k \exp(2\pi i j k / n).\]

The factor of 1/n makes this a true inverse.

In general there are two possible choices for the sign of the exponential in the transform/ inverse-transform pair. GSL follows the same convention as fftpack, using a negative exponential for the forward transform. The advantage of this convention is that the inverse transform recreates the original function with simple Fourier synthesis. Numerical Recipes uses the opposite convention, a positive exponential in the forward transform.

LGSL interface

LGSL provides a simple interface to perform Fourier transforms of real data with the functions fft.fft() and fft.fftinv(). The first function performs the Fourier transform of a column matrix and the second is the inverse Fourier transform.

The function fft.fft() returns a half-complex array. The latter is similar to a column matrix of complex numbers, but it is actually a different object because the numbers are packed together following some specific rules related to the algorithm.

The idea is that you can access the elements of this vector for reading or writing simply by indexing it. You can also obtain the size of the vector using the operator ‘#’. The valid indices for a half-complex object range from 0 to N-1 where N is the size if the vector. Each element of the vector corresponds to the coefficient \(z_k\) defined above.

When performing Fourier transforms, it is important to know that the computation speed can be greatly influenced by the size of the vector. If the size is a power of two, a very efficient algorithm can be used and we can talk in this case of a Fast Fourier Transform (FFT). In addition, the algorithm has the advantage that it does not require any additional workspace. When the size of the vector is not a power of two, we can have two different cases:

• the size is a product of small prime numbers
• the size contains a big (> 7) prime number in its factorization

This detail is important because if the size is a product of small prime numbers, a fast algorithm is still available but it is still somewhat slower and it does require some additional workspace. In the worst case when the size cannot be factorized to small prime numbers, the Fourier transform can still be computed but the calculation is slower, especially for large arrays.

LGSL hides all the details and takes care of choosing the appropriate algorithm based on the size of the vector. It also transparently provides any additional workspace that may be needed for the algorithm. In order to avoid repeated allocation of workspace memory, the workspace allocated is kept in memory and reused if the size of the array does not change. This means that the approach of LGSL is quite optimal if you perform many Fourier transforms (direct or inverse) of the same size.

Even though LGSL takes care of the details automatically, you should be aware of these performance notices because it can make a big difference in real applications. From a practical point of view, it is useful in most cases to always provide samples whose size is a power of two.

Another property of the functions fft.fft() and fft.fftinv() is that they can optionally perform the transformation in place by modifying the original data instead of creating a copy. When a transformation in place is requested, the routine still returns a new vector (either a real matrix or a half-complex array) but this latter will point to the same underlying data of the original vector. The transformation in place can be useful in some cases to avoid unnecessary data copying and memory allocation.

Fourier Transform of Real Data

For real data, the Fourier coefficients satisfy the relation

\[z_k = z_{N-k}^*\]

where N is the size of the vector and k is any integer number from 0 to N-1. Because of this relation, the data is packed in a special type of object called a half-complex array.

To access an element in a half-complex array, you can index it with an integer number between 0 and N-1, inclusive. So, for example:

-- Load all modules in a table called lgsl
lgsl = require("lgsl")

-- get a random number generator
r = lgsl.rng.new()

-- create a vector with random numbers
x = lgsl.matrix.new(256, 1, function() return lgsl.rnd.gaussian(r, 1) end)

-- take the Fourier transform
ft = lgsl.fft.fft(x)

-- print all the coefficients of the Fourier transform
for k=0, #ft-1 do print(ft[k]) end

As shown in the example above, you can use the Lua operator ‘#’ to obtain the size of a half-complex array.

fft.fft(v[, in_place])

Perform the Fourier transform of the real-valued column matrix v. If in_place is true then the original data is altered and the resulting vector will point to the same underlying data of the original vector.

Please note that the value you obtain is not an ordinary matrix but a half-complex array. You can access the elements of such an array by indexing the vector. If you want to have an ordinary matrix you can easily build it with the following instructions:

-- we suppose that ft is an half-complex array
m = lgsl.matrix.cnew(#ft, 1, function(i,j) return ft[i-1] end)

fft.fftinv(hc[, in_place])

Return a column matrix that contains the inverse Fourier transform of the half-complex vector hc. If in_place is true then the original data is altered and the resulting vector will point to the same underlying data of the original vector.

This transformation is the inverse of the function fft.fft(), so that if you perform the two transformations consecutively you will obtain a vector identical to the initial one.

A typical usage of fft.fft_inv() is to revert the transformation made with fft.fft() but by doing some transformations along the way, such as applying a filter. So a typical usage path could be:

-- we assume v is a column matrix with our data
ft = fft.fft(v) -- Fourier transform

-- here we can manipulate the half-complex array 'ft'
-- using the methods `get' and `set'

vt = fft.fftinv(ft) -- we perform the inverse Fourier transform
-- now vt is a vector of the same size of v

FFT example

In this example we will treat a square pulse in the temporal domain. To illustrate a typical example of FFT usage we perform the Fourier Transform of the signal and we cut the higher order frequencies. Then we perform the inverse transform and we compare the result with the original time signal.

So, first we define our square pulse in the time domain. Actually it will be a matrix with just one column:

matrix = require("lgsl.matrix")
fft = require("lgsl.fft")

n, ncut = 256, 16

-- we create a pulse signal in the time domain
y = matrix.new(n, 1, function(i) return  i < n/3 and 0 or (i < 2*n/3 and 1 or 0) end )

If we have graph-toolkit available, we can create two new plots, one for the Fourier transform and one for the signal itself:

graph = require("graph")
pt = graph.plot('Original signal / reconstructed')

pt:addline(graph.filine(function(i) return y[i] end, 1, n), 'black')

Now we are ready to perform:

• the Fourier transform
• cut the higher frequencies
• transform back the signal in the time domain

and plot the results:

ft = fft.fft(y)

complex = require("lgsl.complex")
pf = graph.fibars( function(k) return complex.abs(ft[k]) end, 0, 60)
pf.title = 'FFT Power Spectrum'

for k=ncut, n/2 do ft[k] = 0 end
ytr = fft.fftinv(ft)

pt:addline(graph.filine(function(i) return ytr[i] end, n), 'red')

pt:show()

Fourier transform spectrum

Time signal before (black) and after (red) the transformation

You can observe in the reconstructed signal (the red curve) that we obtain approximately the square pulse, but with a lot of oscillations. These oscillations are an artifact of our transformations, sometimes referred to as the Gibbs phenomenon. The reason is that in order to perfectly reproduce a sharp signal, we would need an infinite number of frequencies in the Fourier transform.