FFTの擬似コード

DFTの回転因子をまとめただけのFFTの擬似コードです。擬似コードというかPythonなので一応実行は出来ますが、きちんとしたFFTの設計は大浦さんのページに詳しく載っています（文字化け注意）。

- [FFT (高速フーリエ・コサイン・サイン変換) の概略と設計法](http://www.kurims.kyoto-u.ac.jp/~ooura/fftman/index.html) 

# 周波数間引き Radix-2 FFT

```math
\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}
```

```py:fft2f.py
import numpy as np
def fft(x):
    n = len(x)
    if n == 1:
        return x
    else:
        # Radix-2
        w  = np.exp(-2j*np.pi/n*np.arange(n//2))
        x0 = x[:n//2] + x[n//2:]
        x1 = x[:n//2] - x[n//2:]
        x0_ = fft( x0 )
        x1_ = fft( x1*w )
        return np.c_[x0_, x1_].ravel()
```

# 時間間引き Radix-2 FFT

```math
\begin{eqnarray}
A_k &=& \sum_{j=0}^{N/2-1} a_{2j} W_{N/2}^{jk} + W^k_N\sum_{j=0}^{N/2-1} a_{2j+1} W_{N/2}^{jk} \\
A_{N/2+k} &=&  \sum_{j=0}^{N/2-1} a_{2j} W_{N/2}^{jk} - W^k_N\sum_{j=0}^{N/2-1} a_{2j+1} W_{N/2}^{jk} 
\end{eqnarray}
```

```py:fft2t.py
import numpy as np
def fft(x):
    n = len(x)
    if n == 1:
        return x
    else:
        # Radix-2
        w  = np.exp(-2j*np.pi/n*np.arange(n//2))
        x0 = fft(x[0::2])
        x1 = fft(x[1::2])*w
        x0_ = x0 + x1
        x1_ = x0 - x1
        return np.r_[x0_, x1_]
```

# 周波数間引き Radix-4 FFT

```math
\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}
```

```py:fft4f.py
import numpy as np
def fft(x):
    n = len(x)
    if n == 1:
        return x
    elif n == 2:
        # Radix-2
        return np.r_[x[0]+x[1], x[0]-x[1]]
    else:
        # Radix-4
        w  = np.exp(-2j*np.pi/n*np.arange(n//4))
        w2 = np.exp(-2j*np.pi/n*2*np.arange(n//4))
        w3 = np.exp(-2j*np.pi/n*3*np.arange(n//4))
        x0 = x[:n//4] + x[n//4*2:n//4*3]
        x1 = x[:n//4] - x[n//4*2:n//4*3]
        x2 = x[n//4:n//4*2] + x[n//4*3:]
        x3 = x[n//4:n//4*2] - x[n//4*3:]
        x0_ = fft( x0 + x2 )
        x1_ = fft( (x1 - 1j*x3)*w )
        x2_ = fft( (x0 - x2)*w2 )
        x3_ = fft( (x1 + 1j*x3)*w3 )
        return np.c_[x0_,
                     x1_,
                     x2_,
                     x3_].ravel()
```

# 時間間引き Radix-4 FFT

```math
\begin{eqnarray}
A_k &=& \sum_{j=0}^{N/4-1} a_{4j} W_{N/4}^{jk} + W^k_N\sum_{j=0}^{N/4-1} a_{4j+1} W_{N/4}^{jk} 
+ W^{2k}_N\sum_{j=0}^{N/4-1} a_{4j+2} W_{N/4}^{jk} + W^{3k}_N\sum_{j=0}^{N/4-1} a_{4j+3} W_{N/4}^{jk} \\
A_{N/4+k} &=& \sum_{j=0}^{N/4-1} a_{4j} W_{N/4}^{jk} - W^{2k}_N\sum_{j=0}^{N/4-1} a_{4j+2} W_{N/4}^{jk}
- i \bigg( W^k_N\sum_{j=0}^{N/4-1} a_{4j+1} W_{N/4}^{jk} - W^{3k}_N\sum_{j=0}^{N/4-1} a_{4j+3} W_{N/4}^{jk} \bigg) \\
A_{2N/4+k} &=& \sum_{j=0}^{N/4-1} a_{4j} W_{N/4}^{jk} - W^k_N\sum_{j=0}^{N/4-1} a_{4j+1} W_{N/4}^{jk} 
+ W^{2k}_N\sum_{j=0}^{N/4-1} a_{4j+2} W_{N/4}^{jk} - W^{3k}_N\sum_{j=0}^{N/4-1} a_{4j+3} W_{N/4}^{jk} \\
A_{3N/4+k} &=& \sum_{j=0}^{N/4-1} a_{4j} W_{N/4}^{jk} - W^{2k}_N\sum_{j=0}^{N/4-1} a_{4j+2} W_{N/4}^{jk}
+ i \bigg( W^k_N\sum_{j=0}^{N/4-1} a_{4j+1} W_{N/4}^{jk} - W^{3k}_N\sum_{j=0}^{N/4-1} a_{4j+3} W_{N/4}^{jk} \bigg)
\end{eqnarray}
```

```py:fft4t.py
import numpy as np
def fft(x):
    n = len(x)
    if n == 1:
        return x
    elif n == 2:
        # Radix-2
        return np.r_[x[0]+x[1], x[0]-x[1]]
    else:
        # Radix-4
        w  = np.exp(-2j*np.pi/n*np.arange(n//4))
        w2 = np.exp(-2j*np.pi/n*2*np.arange(n//4))
        w3 = np.exp(-2j*np.pi/n*3*np.arange(n//4))
        x0 = fft(x[0::4])
        x1 = fft(x[1::4])*w
        x2 = fft(x[2::4])*w2
        x3 = fft(x[3::4])*w3
        x0_ = x0 + x2
        x1_ = x1 + x3
        x2_ = x0 - x2
        x3_ = x1 - x3
        return np.r_[x0_ + x1_,
                     x2_ - 1j*x3_,
                     x0_ - x1_,
                     x2_ + 1j*x3_]
```

# 周波数間引き Split-Radix FFT

```math
\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}
```

```py:fftsf.py
import numpy as np
def fft(x):
    n = len(x)
    if n == 1:
        return x
    elif n == 2:
        # Radix-2
        return np.r_[x[0]+x[1], x[0]-x[1]]
    else:
        # Split-Radix
        w  = np.exp(-2j*np.pi/n*np.arange(n//4))
        w3 = np.exp(-2j*np.pi/n*3*np.arange(n//4))
        x1 = x[:n//4] - x[n//4*2:n//4*3]
        x3 = x[n//4:n//4*2] - x[n//4*3:]
        x0 = fft(x[:n//2] + x[n//2:])
        x1_ = fft( (x1 - 1j*x3)*w )
        x3_ = fft( (x1 + 1j*x3)*w3 )
        return np.c_[x0[0::2],
                     x1_,
                     x0[1::2],
                     x3_].ravel()
```

# 時間間引き Split-Radix FFT

```math
\begin{eqnarray}
A_k &=& \sum_{j=0}^{N/2-1} a_{2j} W_{N/2}^{jk}
+ W^k_N\sum_{j=0}^{N/4-1} a_{4j+1} W_{N/4}^{jk} + W^{3k}_N\sum_{j=0}^{N/4-1} a_{4j+3} W_{N/4}^{jk} \\
A_{N/4+k} &=& \sum_{j=0}^{N/2-1} a_{2j} W_{N/2}^{j(k+N/4)} 
- i \bigg( W^k_N\sum_{j=0}^{N/4-1} a_{4j+1} W_{N/4}^{jk} - W^{3k}_N\sum_{j=0}^{N/4-1} a_{4j+3} W_{N/4}^{jk} \bigg) \\
A_{2N/4+k} &=& \sum_{j=0}^{N/2-1} a_{2j} W_{N/2}^{jk}
- W^k_N\sum_{j=0}^{N/4-1} a_{4j+1} W_{N/4}^{jk} - W^{3k}_N\sum_{j=0}^{N/4-1} a_{4j+3} W_{N/4}^{jk} \\
A_{3N/4+k} &=& \sum_{j=0}^{N/2-1} a_{2j} W_{N/2}^{j(k+N/4)}
+ i \bigg( W^k_N\sum_{j=0}^{N/4-1} a_{4j+1} W_{N/4}^{jk} - W^{3k}_N\sum_{j=0}^{N/4-1} a_{4j+3} W_{N/4}^{jk} \bigg)
\end{eqnarray}
```

```py:fftst.py
import numpy as np
def fft(x):
    n = len(x)
    if n == 1:
        return x
    elif n == 2:
        # Radix-2
        return np.r_[x[0]+x[1], x[0]-x[1]]
    else:
        # Split-Radix
        w  = np.exp(-2j*np.pi/n*np.arange(n//4))
        w3 = np.exp(-2j*np.pi/n*3*np.arange(n//4))
        x0 = fft(x[0::2])
        x1 = fft(x[1::4])*w
        x3 = fft(x[3::4])*w3
        x1_ = x1 + x3
        x3_ = x1 - x3
        return np.r_[x0[:n//4] + x1_,
                     x0[n//4:] - 1j*x3_,
                     x0[:n//4] - x1_,
                     x0[n//4:] + 1j*x3_]
```