Math
As we already know, we can write the Fourier transform in a discrete form like this:
Note, that when working with DFT we use different notation convention. If in a continuous domain we wrote the Fourier transform as \(\hat{f}\), here we notate the sampled (tabbed) function as \(x\) as it’s frequency domain representation as \(X_k\) where \(k\) is the frequency argument.
So, let’s split that heap of sampled function values into two equally sized parts. One part will contain even indexed samples and the other part odd indexed samples, respectively. Then we can rewrite the DFT in the following form:
Now, we want both of this parts to represent a separate self-sufficient DFT transform themselves. To achieve this, we are moving multiplication by 2 from the numerator to denominator, so there it effectively becomes \(N/2\). This way, sigma expression summing from 0 to \(N/2\) becomes a DFT. In the odd part we just moving static stuff outside the sigma operator:
For further convenience we will define left (even) DFT and right (odd) DFT as \(E_k\) and \(O_k\). In this case we can interpret the statement above as:
Now, here is the moment where magic happens: complex exponential functions are periodic, and we can take advantage of that, by reusing already calculated \(E_k\) and \(O_k\) to also calculate \(X_{k + N/2}\). Let’s see, how:
When simplifying this expression, we express exponential degree of sum as a multiplication of two exponential degrees. So, in the first and second sum it just eliminates, because \(e^{-2\pi i m} = 1\) for any \(m\), and the one coefficient on the right will just change it’s sign, because \(e^{-\pi i} = -1\). In detail this would be like:
Which simplifies to:
So, many classical FFT algorithms perform those calculations recursively, reusing the previously computed \(E_k\) and \(O_k\) (which is illustrated by “butterfly” data flow diagram of those algorithms). The complexity of this algorithm family is about O(N log N)