This is a pseudo code of FFT that just summarizes the rotation factors of DFT. Since it is a pseudo code or Python, it can be executed for the time being, but the proper FFT design is described in detail on Mr. Oura's page (be careful of garbled characters).
-Outline and design method of FFT (Fast Fourier Cosine Sine Transform)
Radix-2 FFT
\begin{eqnarray}
A_{2k} &=& \sum_{j=0}^{N/2-1} (a_j + a_{N/2+j}) W_{N/2}^{jk} \\
A_{2k+1} &=& \sum_{j=0}^{N/2-1} ( a_j - a_{N/2+j}) W_N^j W_{N/2}^{jk}
\end{eqnarray}
radix_2_fft.py
from cmath import exp, pi
from itertools import chain
def flatten(l):
return list(chain.from_iterable(l))
def fft(x):
n = len(x)
if(n == 1):
return x
else:
# Radix-2
t1 = (exp(-2j*pi*k/n) for k in xrange(n/2))
x0 = fft( [(a+b) for a, b in zip(x, x[n/2:])] )
x1 = fft( [(a-b)*w for a, b, w in zip(x, x[n/2:], t1)] )
return flatten(zip(x0, x1))
Mixed-Radix FFT (Radix-4 & Radix-2)
\begin{eqnarray}
A_{4k} &=& \sum_{j=0}^{N/4-1} (a_j + a_{N/4+j} + a_{2N/4+j} + a_{3N/4+j}) W_{N/4}^{jk} \\
A_{4k+1} &=& \sum_{j=0}^{N/4-1} (a_j - a_{2N/4+j} - i(a_{N/4+j} - a_{3N/4+j}) ) W_N^j W_{N/4}^{jk} \\
A_{4k+2} &=& \sum_{j=0}^{N/4-1} (a_j - a_{N/4+j} + a_{2N/4+j} - a_{3N/4+j}) W_N^{2j} W_{N/4}^{jk} \\
A_{4k+3} &=& \sum_{j=0}^{N/4-1} (a_j - a_{2N/4+j} + i(a_{N/4+j} - a_{3N/4+j}) ) W_N^{3j} W_{N/4}^{jk}
\end{eqnarray}
mixed_radix_fft.py
from cmath import exp, pi
from itertools import chain
def flatten(l):
return list(chain.from_iterable(l))
def fft(x):
n = len(x)
if(n == 1):
return x
elif(n == 2):
# Radix-2
return [x[0]+x[1], x[0]-x[1]]
else:
# Radix-4
t1 = (exp(-2j*pi*k/n) for k in xrange(n/4))
t2 = (exp(-2j*pi*2*k/n) for k in xrange(n/4))
t3 = (exp(-2j*pi*3*k/n) for k in xrange(n/4))
x0 = fft( [(a+b+c+d) for a, b, c, d in zip(x, x[n/4:], x[n*2/4:], x[n*3/4:])] )
x1 = fft( [(a-c-1j*(b-d))*w for a, b, c, d, w in zip(x, x[n/4:], x[n*2/4:], x[n*3/4:], t1)] )
x2 = fft( [(a-b+c-d)*w for a, b, c, d, w in zip(x, x[n/4:], x[n*2/4:], x[n*3/4:], t2)] )
x3 = fft( [(a-c+1j*(b-d))*w for a, b, c, d, w in zip(x, x[n/4:], x[n*2/4:], x[n*3/4:], t3)] )
return flatten(zip(x0, x1, x2, x3))
Split-Radix FFT
\begin{eqnarray}
A_{2k} &=& \sum_{j=0}^{N/2-1} (a_j + a_{N/2+j}) W_{N/2}^{jk} \\
A_{4k+1} &=& \sum_{j=0}^{N/4-1} (a_j - a_{2N/4+j} - i(a_{N/4+j} - a_{3N/4+j}) ) W_N^j W_{N/4}^{jk} \\
A_{4k+3} &=& \sum_{j=0}^{N/4-1} (a_j - a_{2N/4+j} + i(a_{N/4+j} - a_{3N/4+j}) ) W_N^{3j} W_{N/4}^{jk}
\end{eqnarray}
split_radix_fft.py
from cmath import exp, pi
from itertools import chain
def flatten(l):
return list(chain.from_iterable(l))
def fft(x):
n = len(x)
if(n == 1):
return x
elif(n == 2):
# Radix-2
return [x[0]+x[1], x[0]-x[1]]
else:
# Split-Radix
t1 = (exp(-2j*pi*k/n) for k in xrange(n/4))
t2 = (exp(-2j*pi*3*k/n) for k in xrange(n/4))
x0 = fft( [(a+b) for a, b in zip(x, x[n/2:])] )
x1 = fft( [(a-c-1j*(b-d))*w for a, b, c, d, w in zip(x, x[n/4:], x[n*2/4:], x[n*3/4:], t1)] )
x2 = fft( [(a-c+1j*(b-d))*w for a, b, c, d, w in zip(x, x[n/4:], x[n*2/4:], x[n*3/4:], t2)] )
return flatten(zip(x0[0::2], x1, x0[1::2], x2))