Discrete Fourier Transform (DFT) takes $x_0,\dots,x_{N-1}$ and produces $X_0,\dots,X_{N-1}$ by the following equality

$$X_k = \sum^{N-1}_{n=0} x_n e^{-i2\pi kn/N}.$$

Fast Fourier Transform (FFT) does this computation effectively in $O(N \log N)$ time. Possibility of doing the computation in a more effective way raises from decomposition of the values into odd and even elements, and applying the computation recursively.

Let $\omega_N=e^{-2i\pi /N}$, so we have

$$ X_k = \sum^{N-1}{n=0} x_n \omega^{nk}N = \sum^{N/2-1}{n=0} x{2n} \omega^{2nk}_N

  • \sum^{N/2-1}{n=0} x{2n+1} \omega^{(2n+1)k}N = $$ $$ = \sum^{N/2-1}{n=0} x_{2n} \omega^{2nk}_N
  • \omega^k_N \sum^{N/2-1}{n=0} x{2n+1} \omega^{2nk}_N. $$

We note that $\omega^{2nN/2}N=e^{-2i\pi(2nN/2)/N}=e^{-2i\pi n}=1$ (because of $e^{i\pi}=\cos(\pi)+i\sin(\pi)=-1$ by definition). The element $$ X{k+N/2} = \sum^{N/2-1}{n=0} x{2n} \omega^{2n(k+N/2)}_N

  • \omega^{k+N/2}N \sum^{N/2-1}{n=0} x_{2n+1} \omega^{2n(k+N/2)}N = $$ and by a nice sequence of simplifications we have $$ = \sum^{N/2-1}{n=0} x_{2n} \omega^{2nk}_N
  • (-1) \omega^k_N \sum^{N/2-1}{n=0} x{2n+1} \omega^{2nk}_N, $$ which the same equation as for $X_k$ up to the $(-1)$.

Code

Disclaimer: this code was provided by Morass for the acm-solutions for BI-ACM course on FIT CTU.

This code is made for $|{\rm src}|=2^i$ and will probably not work on source arrays of other sizes.

typedef complex<double> cpx;

void fft(const cpx *src, cpx *res, int n, const cpx &wn, cpx *tmp){
    if(n==1){ *res = *src; return; }
    int N = n/2;
    cpx *osrc=tmp, *esrc=tmp+N, *ores=tmp+2*N, *eres=tmp+3*N;
    for(int i=0; i<N; ++i){
        esrc[i] = src[2*i];
        osrc[i] = src[2*i+1];
    }
    fft(esrc, eres, N, (wn*wn), tmp+4*N);
    fft(osrc, ores, N, (wn*wn), tmp+4*N);
    cpx w(1, 0);
    for(int i=0; i<N; ++i){
        res[i  ] = eres[i]+w*ores[i];
        res[i+N] = eres[i]-w*ores[i];
        w *= wn;
    }
}

The implementation uses the class in the c++ std #include <complex>.

typedef complex<double> cpx;

In case whre the array size is one, the solution is trivial.

    if(n==1){ *res = *src; return; }

Otherwise, we split the problem in-two: odd and even-indexed subproblems. The size of these sub-problems is half of the original.

    int N = n/2;

The computation for the subproblems will happen in a pre-allocated temporary array. Each recursive call has problem of half the size and the same-level subproblems reuse the memory. We start with empty tmp, and require 4*|src|/2 for first subproblems, which implies that 4*|src| size suffices for the tmp array.

    cpx *osrc=tmp, *esrc=tmp+N, *ores=tmp+2*N, *eres=tmp+3*N;

We prepare the input for the subproblem by splitting our input into even and odd-indexed values.

    for(int i=0; i<N; ++i){
        esrc[i] = src[2*i];
        osrc[i] = src[2*i+1];
    }

And call the recursion, note how this call relates to the $\sum^{N/2-1}{n=0} x{2n} \omega^{2nk}N$ and $\sum^{N/2-1}{n=0} x_{2n+1} \omega^{2nk}_N$ respectively.

  • We pass arrays with relevant input and output, which are half of the original’s size,
  • the passed N is n/2, indicating size of the problem is halved,
  • by passing (wn*wn) we set $\omega:=\omega^2$, as seen in the equation,
  • last, tmp+4N is the first non-used element of the tmp array, and is always available by the argument above.
    fft(esrc, eres, N, (wn*wn), tmp+4*N);
    fft(osrc, ores, N, (wn*wn), tmp+4*N);

The solution needs to be put together from the sub-solutions according to the equations for $X_k$ and $X_{k+N/2}$.

        res[i  ] = eres[i]+w*ores[i];
        res[i+N] = eres[i]-w*ores[i];
        w *= wn;