What does Cooley–Tukey FFT algorithm calculates in order to compute DFT
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 &amp;=&amp; \sum_{n = 0}^{N-1} x_n e^{-2j\pi k n/N} \\
&amp;=&amp; \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 }} \\
 &amp;=&amp; \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 &lt;complex&gt;
#include &lt;iostream&gt;
#include &lt;valarray&gt;
 
const double PI = 3.141592653589793238460;
 
typedef std::complex&lt;double&gt; Complex;
typedef std::valarray&lt;Complex&gt; CArray;
 
// Cooley–Tukey FFT (in-place, divide-and-conquer)
// Higher memory requirements and redundancy although more intuitive
void fft(CArray&amp; x)
{
    const size_t N = x.size();
    if (N &lt;= 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 &lt; 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;
    }
}