| tags:[ programming algorithm ]

# Fast Fourier Transform (FFT) c++ code explained

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;
```