DFTの回転因子をまとめただけのFFTの擬似コードです。擬似コードというかPythonなので一応実行は出来ますが、きちんとしたFFTの設計は大浦さんのページに詳しく載っています（文字化け注意）。
周波数間引き 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}
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
\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}
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
\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}
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
\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}
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
\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}
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
\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}
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_]