[PYTHON] FFT pseudo code

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))