What does Cooley–Tukey FFT algorithm calculates in order to compute DFT

0

I'm relative new to this subject, I've watched many videos explaining FFT and DFT and read some articles.

I wanted to see how I could implement FFT in C++ and then I found this code, it works but I don't fully understand it, for example, I'm not sure what is he using the nested-for loops for, I thought it's for Matrix multiplication but maybe it's just using the DFT formula for each element of the series, in addition why does he needs to reverse the bits?

if someone could explain in "simple words" the idea behind this code I'll appriciate it because I just started learning it and I've come crossed super complicated stuff which I don't understand.

Here is the FFT Function I'm using:

void fft(Iter_T a, Iter_T b, int log2n)
{
typedef typename iterator_traits<Iter_T>::value_type complex;
const complex J(0, 1);
int n = 1 << log2n;
for (unsigned int i = 0; i < n; ++i) {
b[ReverseBits(i, log2n)] = a[i];
}
for (int s = 1; s <= log2n; ++s) {
int m = 1 << s;
int m2 = m >> 1;
complex w(1, 0);
complex wm = exp(-J * (PI / m2));
for (int j = 0; j < m2; ++j) {
for (int k = j; k < n; k += m) {
complex t = w * b[k + m2];
complex u = b[k];
b[k] = u + t;
b[k + m2] = u - t;
}
w *= wm;
}
}
}
• There are tons of websites and videos that explain this. What have you read and what do you not understand ? May 12 '21 at 13:43

This is the iterative version of the algorithm. Most of the difficulty there is tracing the indices of the subarrays in each stage of the algorithm. For instance bit reversing is used so that the even inputs are moved to the first half of the subarray and odd inputs are moved to the end of the subarray at all levels, to improve memory access locality.

Conceptually the Cooley-Tukey is a recursive algorithm in which each step you split the input in even/odd indices subarrays, and compute the first and second half of the DFT.

$$\begin{eqnarray} X_k &=& \sum_{n = 0}^{N-1} x_n e^{-2j\pi k n/N} \\ &=& \underbrace{\sum_{n = 0}^{N/2-1} x_{2n} e^{-2j\pi k\, 2n/N}}_{\textrm{even term}} + \underbrace{\sum_{n = 0}^{N-1} x_{2n+1} e^{-2j\pi k\, (2n+1)/N}}_{\textrm{odd term }} \\ &=& \underbrace{\sum_{n = 0}^{N/2-1} x_{2n} e^{-2j\pi k\, n/(N/2)}}_{E_k} + e^{2j\pi\,k/N}\underbrace{\sum_{n = 0}^{N-1} x_{2n+1} e^{-2j\pi k\, n/(N/2)}}_{O_k} \\ \end{eqnarray}$$

Then $$E_k$$ and $$O_k$$ are themselves the coefficients of a DFT of the decimated even/odd indexed terms of $$x$$ (of size N/2).

The loop computes two terms one for $$X_k = E_k + O_k$$ and, (using the periodicity $$E_k = E_{k+N/2}$$ and $$O_k=O_{k+N/2}$$), $$X_{k+N/2} = E_k + (-1)O_k$$ (the $$-1$$ comes from the exponential mulitplying the $$O_k$$).

Here is a more readable C++ implementation (from rosetta code)

#include <complex>
#include <iostream>
#include <valarray>

const double PI = 3.141592653589793238460;

typedef std::complex<double> Complex;
typedef std::valarray<Complex> CArray;

// Cooley–Tukey FFT (in-place, divide-and-conquer)
// Higher memory requirements and redundancy although more intuitive
void fft(CArray& x)
{
const size_t N = x.size();
if (N <= 1) return;

// divide
CArray even = x[std::slice(0, N/2, 2)];
CArray  odd = x[std::slice(1, N/2, 2)];

// conquer
fft(even);
fft(odd);

// combine
for (size_t k = 0; k < N/2; ++k)
{
Complex t = std::polar(1.0, -2 * PI * k / N) * odd[k];
x[k    ] = even[k] + t;
x[k+N/2] = even[k] - t;
}
}
• I tried using the rosetta implementation with the inputs from the code I'm using at the moment and it computes different outputs, any idea why? I am using these inputs: const Complex test[] = { 0.0, 1.1, 3.3, 4.4, 4.4, 3.3, 1.1, 0.0 }; May 12 '21 at 14:15
• Check the order of the coefficients, [17.6, -9.0669 + 3.7556i, 0, 0.2669 - 0.6444i, 0, ...], the the other three are symmetric, complex conjugates, or normalization factors. Also you could have a reference implementation that just computes the coefficients by the definition with O(N^2) complexity to help you debugging.
– Bob
May 12 '21 at 16:14