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;