OpenCL Fast Fourier Transform

Eric Bainville - May 2010

Fast Fourier Transform algorithms


The Fast Fourier Transform is a family of algorithms computing the Discrete Fourier Transform (DFT) in time O(N.log(N)). The most commonly used FFT algorithms use a divide-and-conquer approach similar to the algorithm of Cooley and Tukey (1965). The Wikipedia FFT page is a good starting point for the interested reader.

The computation of a DFT of length N is done by splitting the input sequence into a fixed small number of smaller subsequences, compute their DFT, and assemble the outputs to build the final sequence. Providing the split and assembly steps have linear cost, the complexity of the algorithm will be O(N.log(N)).

One of the simplest subdivision is to take the sub-sequence xe with even indices 0,2,4,...,N-2, and xo with odd indices 1,3,...,N-1, supposing N is even. Compute the DFT of each sub-sequence, ye=DFTN/22,xe) and yo=DFTN/22,xo). Then we can build y=DFTN(θ,x) as follows:
yk = yek + θk yok, and
yk+N/2 = yek - θk yok, for k=0..N/2-1.

More details on subdivisions

In this section, we enter into more details about the various subdivision schemes used in FFT, known as decimation in frequency (DIF), and decimation in time (DIT). Both actually refer to the same algorithm.

General case of FFT subdivision.

Consider a sequence x=(x0,x1,...,xN-1) of length N, and a factorization N=P*Q. We denote the DFT of x of length N by DFTN(x) or DFTN(x0,x1,...,xN-1). Let θ=e-2*i*Π/N, y=DFTN(x) is by definition:

yk = θ0.k.x0 + θ1.k.x1 + θ2.k.x2 + ... (N terms).

We consider the P subsequences of x with length Q defined by elements with the same index modulo P:
x[u]=(xu,xu+P,xu+2.P,...,xu+(Q-1).P) with u=0..P-1.

Highlighting the terms of x[u] in yk, we get:
yk = ... + θu.k.xu + ... + θ(u+P).k.xu+P + ... + θ(u+2.P).k.xu+2.P + ..., or
yk = Σu θu.k (xu + θP.k.xu+P + θ2.P.k.xu+2.P + ...), and finally
yk = Σu θu.k DFTQ(x[u])k, where element k of a DFT of length Q is taken modulo Q.

Let's group together the P values of k using the same element v of y[u]=DFTQ(x[u]), by expressing k as k=v+Q.j, with v=0..Q-1 and j=0..P-1:
yv+Q.j = Σu θu.(v+Q.j) y[u]v, or
yv+Q.j = Σu θQ.u.ju.v y[u]v), and finally
yv+Q.j = DFTP0.v y[0]v1.v y[1]v,...,θ(P-1).v y[P-1]v)j.

We have shown that the DFTN can be computed by:
- compute P×DFTQ,
- scale all the N resulting values by the "twiddle factors" θu.v,
- compute Q×DFTP.

Decimation in time (DIT, left) and decimation in frequency (DIF, right). In blue the input sequence, and in dark blue one of the P sequences of length Q input of DFTQ. In orange the output sequence, and in dark orange one of the Q sequences of length P output of DFTP.

This algorithm is called either radix-P decimation in time (DIT), or radix-Q decimation in frequency (DIF). Usually, the name is chosen according to the smallest value between P and Q. A radix-2 DIT algorithm corresponds to P=2 and Q=N/2; it starts with 2 DFTN/2. A radix-2 DIF algorithm corresponds to P=N/2 and Q=2; it starts with N/2 DFT2.

In the next page, Reference implementations, we show how the performance of the reference implementations is measured.