「FFT ってどんなんだったっけ？」となったときに思い出すための備忘録。
偶関数と奇関数
偶関数は
f(x) = f(-x)
となる関数。
\cos{x} = \cos{\left(-x\right)}
奇関数は
\displaylines {
-f(x) = f(-x) \\
\int_{-L}^L f \left( \theta \right) \ d \theta = 0
}
となる関数。
- \sin{x} = \sin{\left(-x\right)}
偶関数と奇関数の積は奇関数になるので
\displaylines {
\int_{-\pi}^{\pi} \cos{m \theta} \sin{n \theta} \ d \theta = 0 \\
\left( m,n は整数とする \right)
}
偶関数と偶関数の積は偶関数で
\displaylines {
\frac{1}{\pi} \int_{-\pi}^{\pi} { \cos{m \theta} \cos{n \theta} \ d \theta }
= \frac{1}{\pi} \int_{0}^{2 \pi} \cos{m \theta} \cos{n \theta} \ d \theta = \delta_{m n} \\
\left( m,n は整数で、\delta_{mn} はクロネッカーのデルタ \right)
}
奇関数と奇関数の積は偶関数で
\displaylines {
\frac{1}{\pi} \int_{-\pi}^{\pi} { \sin{m \theta} \sin{n \theta} \ d \theta }
= \frac{1}{\pi} \int_{0}^{2 \pi} { \sin{m \theta} \sin{n \theta} \ d \theta } = \delta_{m n} \\
\left( m,n は整数で、\delta_{mn} はクロネッカーのデルタ \right)
}
になります。
クロネッカーのデルタ $\delta_{mn}$ は、プログラムにすると ((m == n) ? 1 : 0) です。
フーリエ級数
関数を三角関数を使って表現
関数 $f \left( \theta \right)$ が $\left( -\pi \le \theta \le \pi \right)$ の範囲で
f \left( \theta \right) = C + \sum_{n=1}^{\infty} \left( a_n \cos{n \theta} + b_n \sin{n \theta} \right)
にできる周期関数として
\begin{array}{l}
\int_{-\pi}^{\pi} { f \left( \theta \right) \cos{\theta} \ d \theta } \\
\int_{-\pi}^{\pi} { f \left( \theta \right) \sin{\theta} \ d \theta } \\
\int_{-\pi}^{\pi} { f \left( \theta \right) \ d \theta } \\
\end{array}
を求めると
\begin{eqnarray}
\int_{-\pi}^{\pi} { f \left( \theta \right) \cos{\theta} \ d \theta }
& = & \left( C \int_{-\pi}^{\pi} { \cos{\theta} \ d \theta } = 0 \right) \\
& + & \left( \left( a_1 \int_{-\pi}^{\pi} { \cos{\theta} \cos{\theta} \ d \theta } = a_1 \pi \right) \right)
+ \left( a_2 \int_{-\pi}^{\pi} { \cos{2 \theta} \cos{\theta} \ d \theta } = 0 \right) + \cdots \\
& + & \left( b_1 \int_{-\pi}^{\pi} { \sin{\theta} \cos{\theta} \ d \theta } = 0 \right)
+ \left( b_2 \int_{-\pi}^{\pi} { \sin{2 \theta} \cos{\theta} \ d \theta } = 0 \right) + \cdots \\
& = & a_1 \pi \\
\int_{-\pi}^{\pi} { f \left( \theta \right) \sin \theta \ d \theta }
& = & \left( C \int_{-\pi}^{\pi} { \sin{\theta} \ d \theta } = 0 \right) \\
& + & \left( a_1 \int_{-\pi}^{\pi} { \cos{\theta} \sin{\theta} \ d \theta } = 0 \right)
+ \left( a_2 \int_{-\pi}^{\pi} { \cos{2 \theta} \sin{\theta} \ d \theta } = 0 \right) + \cdots \\
& + & \left( \left( b_1 \int_{-\pi}^{\pi} { \sin{\theta} \sin{\theta} \ d \theta } = b_1 \pi \right) \right)
+ \left( b_2 \int_{-\pi}^{\pi} { \sin{2 \theta} \sin{\theta} \ d \theta } = 0 \right) + \cdots \\
& = & b_1 \pi \\
\int_{-\pi}^{\pi} { f \left( \theta \right) \ d \theta }
& = & \left( \left( C \int_{-\pi}^{\pi} { d \theta } = 2 C \pi \right) \right) \\
& + & \left( a_1 \int_{-\pi}^{\pi} { \cos{\theta} \ d \theta } = 0 \right)
+ \left( a_2 \int_{-\pi}^{\pi} { \cos{2 \theta} \ d \theta } = 0 \right)
+ \cdots \\
& + & \left( b_1 \int_{-\pi}^{\pi} { \sin{\theta} \ d \theta } = 0 \right)
+ \left( b_2 \int_{-\pi}^{\pi} { \sin{2 \theta} \ d \theta } = 0 \right)
+ \cdots \\
& = & 2 C \pi \\
\end{eqnarray}
二重括弧$((\cdots))$以外はゼロになります。同様に
\begin{array}{l}
\int_{-\pi}^{\pi} { f \left( \theta \right) \cos{n \theta} \ d \theta } \\
\int_{-\pi}^{\pi} { f \left( \theta \right) \sin{n \theta} \ d \theta } \\
\end{array}
では
\begin{eqnarray}
\int_{-\pi}^{\pi} { f \left( \theta \right) \cos{n \theta} \ d \theta }
& = & \left( C \int_{-\pi}^{\pi} { \cos{n \theta} \ d \theta } = 0 \right) \\
& + & \cdots
+ \left( \left( a_n \int_{-\pi}^{\pi} { \cos{n \theta} \cos{n \theta} \ d \theta } = a_n \pi \right) \right)
+ \cdots \\
& + & \cdots
+ \left( b_m \int_{-\pi}^{\pi} { \sin{m \theta} \cos{n \theta} \ d \theta } = 0 \right)
+ \cdots \\
& = & a_n \pi \\
\int_{-\pi}^{\pi} { f \left( \theta \right) \sin{n \theta} \ d \theta }
& = & \left( C \int_{-\pi}^{\pi} { \sin{n \theta} \ d \theta } = 0 \right) \\
& + & \cdots
+ \left( a_m \int_{-\pi}^{\pi} { \cos{m \theta} \sin{n \theta} \ d \theta } = 0 \right)
+ \cdots \\
& + & \cdots
+ \left( \left( b_n \int_{-\pi}^{\pi} { \sin{n \theta} \sin{n \theta} \ d \theta } = b_n \pi \right) \right)
+ \cdots \\
& = & b_n \pi \\
\end{eqnarray}
ゼロの項を省いて整理すると
\begin{eqnarray}
a_0 \pi = 2 C \pi & = & \int_{-\pi}^{\pi} { f \left( \theta \right) \cos{0} \ d \theta } = \int_{-\pi}^{\pi} { f \left( \theta \right) \ d \theta }\\
b_0 \pi & = & \int_{-\pi}^{\pi} { f \left( \theta \right) \sin{0} \ d \theta } = 0 \\
a_n \pi & = & \int_{-\pi}^{\pi} { f \left( \theta \right) \cos{n} \theta \ d \theta } \\
b_n \pi & = & \int_{-\pi}^{\pi} { f \left( \theta \right) \sin{n} \theta \ d \theta } \\
C & = & \frac{a_0}{2} \\
\end{eqnarray}
だから
\begin{eqnarray}
a_n & = & \frac{1}{\pi} \int_{-\pi}^{\pi} { f \left( \theta \right) \cos{n \theta} \ d \theta } \\
b_n & = & \frac{1}{\pi} \int_{-\pi}^{\pi} { f \left( \theta \right) \sin{n \theta} \ d \theta } \\
f \left( \theta \right) & = & \frac{a_0}{2} + \sum_{n=1}^{\infty} \left( a_n \cos{n \theta} + b_n \sin{n \theta} \right)
\end{eqnarray}
になにります。
三角関数から指数関数(複素数)へ
添字 $n$ を負の方向に拡張して
\begin{eqnarray}
a_{-n} & = & \frac{1}{\pi} \int_{-\pi}^{\pi} { f \left( \theta \right) \cos{\left( -n \theta \right)} \ d \theta } \\
-b_{-n} & = & - \frac{1}{\pi} \int_{-\pi}^{\pi} { f \left( \theta \right) \sin{\left( -n \theta \right)} \ d \theta }
\end{eqnarray}
これを後で使います。
三角関数を指数関数表現
\begin{eqnarray}
\cos{n \theta} & = & \frac{ e^{i\ n \theta} + e^{-i\ n \theta} }{2} \\
\sin{n \theta} & = & \left( -i \right) \frac{ e^{i\ n \theta} - e^{-i\ n \theta} }{2} \\
\end{eqnarray}
にして、$f \left( \theta \right)$ に代入すると
\begin{eqnarray}
f \left( \theta \right) & = & \frac{a_0}{2} + \sum_{n=1}^{\infty} \left( a_n \cos{n \theta} + b_n \sin{n \theta} \right) \\
& = & \frac{a_0}{2} + \sum_{n=1}^{\infty} \left( a_n \frac{ e^{i\ n \theta} + e^{-i\ n \theta} }{2} - i\ b_n \frac{ e^{i\ n \theta} - e^{-i\ n \theta} }{2} \right) \\
& = & \frac{a_0}{2} + \sum_{n=1}^{\infty} \left( \frac{a_n - i\ b_n}{2} e^{i\ n \theta} + \frac{ a_n + i\ b_n }{2} e^{-i\ n \theta} \right)
\end{eqnarray}
ですが、ここで負の添字と $\left(b_0=0\right)$ を使うと
\begin{eqnarray}
f \left( \theta \right) & = & \frac{a_0 - i\ b_0}{2} + \sum_{n=1}^{\infty} \left( \frac{a_n - i\ b_n}{2} e^{i\ n \theta} + \frac{ a_{-n} - i\ b_{-n} }{2} e^{i \left(-n \theta \right)} \right) \\
& = & \sum_{n=-\infty}^{-1} \left( \frac{ a_n - i\ b_n }{2} e^{i\ n \theta } \right) + \left( \frac{a_0 - i\ b_0}{2} e^0 \right) + \sum_{n=1}^{\infty} \left( \frac{a_n - i\ b_n}{2} e^{i\ n \theta} \right) \\
& = & \sum_{n=-\infty}^{\infty} { \frac{a_n - i\ b_n}{2} e^{i\ n \theta} }
\end{eqnarray}
とできます。さらに
c_n = \frac{a_n - i\ b_n}{2}
とすると
\begin{eqnarray}
c_n & = & \frac{a_n - i\ b_n}{2} \\
& = & \frac{1}{2} \left( \frac{1}{\pi} \int_{-\pi}^{\pi} { f \left( \theta \right) \cos{n \theta} \ d\theta} - i \frac{1}{\pi} \int_{-\pi}^{\pi} { f \left( \theta \right) \sin{n \theta} \ d\theta} \right) \\
& = & \frac{1}{2\pi} \int_{-\pi}^{\pi} { f \left( \theta \right) \left( \cos{n \theta} - i \sin{n \theta} \right) \ d\theta} \\
& = & \frac{1}{2\pi} \int_{-\pi}^{\pi} { f \left( \theta \right) \{ \cos{\left( -n \theta \right)} + i \sin{\left( -n \theta \right)} \} \ d\theta} \\
& = & \frac{1}{2\pi} \int_{-\pi}^{\pi} { f \left( \theta \right) e^{-i\ n \theta} \ d\theta}
\end{eqnarray}
だから
\begin{eqnarray}
f \left( \theta \right) & = & \sum_{n=-\infty}^{\infty} { c_n e^{i\ n \theta} } \\
c_n & = & \frac{1}{2\pi} \int_{-\pi}^{\pi} { f \left( \theta \right) e^{-i\ n \theta} \ d\theta}
\end{eqnarray}
が得られます。
周期の変更
関数 $f \left( \theta \right)$ の周期を範囲 $\left( -\pi \le \theta \le \pi \right)$ から $\left( -L \le x \le L \right)$ にするため
\begin{eqnarray}
x & = & \frac{L}{\pi} \theta \\
\theta & = & \frac{\pi}{L} x \\
d \theta & = & \frac{\pi}{L} dx
\end{eqnarray}
により
\begin{eqnarray}
f \left( \theta \right) & = & f \left( \frac{\pi}{L} x \right) = \sum_{n=-\infty}^{\infty} { c_n e^{i \frac{\pi}{L} n x} } = f_L (x) \\
a_n & = & \frac{1}{\pi} \int_{-\pi}^{\pi} { f \left( \theta \right) \cos{n \theta} \ d \theta }
= \frac{1}{L} \int_{-L}^{L} { f_L \left( x \right) \cos{\left( \frac{\pi}{L}n x \right)} \ d x } \\
b_n & = & \frac{1}{\pi} \int_{-\pi}^{\pi} { f \left( \theta \right) \sin{n \theta} \ d \theta }
= \frac{1}{L} \int_{-L}^{L} { f_L \left( x \right) \sin{\left( \frac{\pi}{L}n x \right)} \ d x } \\
c_n & = & \frac{1}{2\pi} \int_{-\pi}^{\pi} { f \left( \theta \right) e^{-i\ n \theta} \ d\theta}
= \frac{1}{2 L} \int_{-L}^{L} { f_L \left( x \right) e^{-i \frac{\pi}{L} n x} \ d x}
\end{eqnarray}
とし、ここから $f(x)$ を
f_L \left( x \right) → f \left( x \right) \cdots \left( -L \le x \le L \right)
に置き換えて
\begin{eqnarray}
f \left( x \right) & = & \sum_{n=-\infty}^{\infty} { c_n e^{i \frac{\pi}{L} n x} } \\
a_n & = & \frac{1}{L} \int_{-L}^{L} { f \left( x \right) \cos{\left( \frac{\pi}{L}n x \right)} \ d x } \\
b_n & = & \frac{1}{L} \int_{-L}^{L} { f \left( x \right) \sin{\left( \frac{\pi}{L}n x \right)} \ d x } \\
c_n & = & \frac{1}{2 L} \int_{-L}^{L} { f \left( x \right) e^{-i \frac{\pi}{L} n x} \ d x} \\
\end{eqnarray}
を使います。
周期の範囲を無限大へ
関数 $f \left( x \right)$ の周期を範囲 $\left( -L \le x \le L \right)$ から $\left( -\infty \le x \le \infty \right)$ にして、周期性のない関数も対象にします。
\begin{eqnarray}
f \left( x \right) & = & \sum_{n=-\infty}^{\infty} { c_n e^{i \frac{\pi}{L} n x} } \\
c_n & = & \frac{1}{2 L} \int_{-L}^{L} { f \left( x \right) e^{-i \frac{\pi}{L} n x} \ d x}
\end{eqnarray}
周期範囲の係数 $L$ を { $L → rL$ } に変更して
\begin{eqnarray}
f \left( x \right) & = & \sum_{n=-\infty}^{\infty} { c_n e^{i \frac{\pi}{r L} n x} } \\
c_n & = & \frac{1}{2 r L} \int_{-r L}^{r L} { f \left( x \right) e^{-i \frac{\pi}{r L} n x} \ d x}
\end{eqnarray}
$r$ が増えても、$f \left( x \right)$ の範囲 $\left( -L \le x \le L \right)$ に変化がなく、それ以外の範囲で $f \left( x \right) = 0$ になる(単パルスの)場合、$r → \infty$ とすると
f \left( x \right) = \int_{-\infty}^{\infty}{ c_n e^{i \frac{\pi}{r L} n x} d n }
になります。$r=1$ の Σ では、幅$1$、高さ $c_n$ の面積の総和ですが、$r$ が増やすと、幅 $1/r$ 倍、高さ $r$ 倍したグラフでは、分解能が高くなっていき、$r → \infty$ では連続したグラフになります。
\begin{eqnarray}
\omega & = & \frac{\pi}{r L} n \\
d \omega & = & \frac{\pi}{r L} d n \\
d n & = & \frac{r L}{\pi} d \omega
\end{eqnarray}
として、$f \left( x \right)$ の $c_n$ を展開すると
\begin{eqnarray}
f \left( x \right) & = & \int_{-\infty}^{\infty}{ \left( \frac{1}{2 r L} \int_{-r L}^{r L} { f \left( x \right) e^{-i \frac{\pi}{r L} n x} \ d x} \right) e^{i \frac{\pi}{r L} n x} d n } \\
& = & \int_{-\infty}^{\infty}{ \left( \frac{1}{2 r L} \int_{-r L}^{r L} { f \left( x \right) e^{-i \omega x} \ d x} \right) e^{i \omega x} \frac{r L}{\pi} d \omega } \\
& = & \int_{-\infty}^{\infty}{ \left( \frac{1}{2 \pi} \int_{-r L}^{r L} { f \left( x \right) e^{-i \omega x} \ d x} \right) e^{i \omega x} d \omega }
\end{eqnarray}
$r → \infty$ として整理すると
f \left( x \right) = \frac{1}{ \sqrt{2 \pi} } \int_{-\infty}^{\infty}{ \left( \frac{1}{ \sqrt{2 \pi} } \int_{-\infty}^{\infty} { f \left( x \right) e^{-i \omega x} \ d x} \right) e^{i \omega x} d \omega }
内側の積分がフーリエ変換
\mathcal{F}[f \left( x \right)] = F \left( \omega \right) = \frac{1}{ \sqrt{2 \pi} } \int_{-\infty}^{\infty} { f \left( x \right) e^{-i \omega x} \ d x}
外側の積分が逆フーリエ変換
\mathcal{F}^{-1}[F \left( \omega \right)] = f \left( x \right) = \frac{1}{ \sqrt{2 \pi} } \int_{-\infty}^{\infty} { F \left( \omega \right) e^{i \omega x} \ d \omega}
になります。内外の違いは、指数関数の指数部の符号だけ。
離散フーリエ変換
サンプリング データをフーリエ変換する。
基本計算
フーリエ級数まで戻って
\begin{eqnarray}
f \left( x \right) & = & \sum_{n=-\infty}^{\infty} { c_n e^{i \frac{\pi}{L} n x} } \\
c_n & = & \frac{1}{2 L} \int_{-L}^{L} { f \left( x \right) e^{-i \frac{\pi}{L} n x} \ d x}
\end{eqnarray}
周期の範囲 $\left( -\pi \le \theta \le \pi \right)$ を $\left( 0 \le \theta \le 2 \pi \right)$ に変更しても同じだから、$\left( -L \le x \le L \right)$ は $\left( 0 \le x \le 2 L \right)$ に変更して、$\lambda = 2 L$ とすると
\begin{eqnarray}
\lambda & = & 2 L \\
f \left( x \right) & = & \sum_{n=-\infty}^{\infty} { c_n e^{i \frac{2 \pi}{\lambda} n x} } \\
c_n & = & \frac{1}{\lambda} \int_{0}^{\lambda} { f \left( x \right) e^{-i \frac{2 \pi}{\lambda} n x} \ d x}
\end{eqnarray}
とできます。データを処理するには、積分を(幅1の棒グラフ)総和にします。サンプリング数を $\left( N = \lambda \right)$ とすると
c_n = \frac{1}{N} \sum_{x=0}^{N-1} { f \left( x \right) e^{-i \frac{2 \pi}{N} n x} }
になります。$c_n$ は、データの分解能 $\left( N = \lambda \right)$ を超えたものを無視して、$\left( 0 \le n \lt N \right)$ とすると、$f \left( x \right)$ は
f \left( x \right) = \sum_{n=0}^{N-1} { c_n e^{i \frac{2 \pi}{N} n x} }
となります。これで、フーリエ変換と逆変換
\begin{eqnarray}
f \left( v \right) & = & \sum_{u=0}^{N-1} { F \left( u \right) e^{i \frac{2 \pi}{N} u v} } \\
F \left( u \right) & = & \frac{1}{N} \sum_{v=0}^{N-1} { f \left( v \right) e^{-i \frac{2 \pi}{N} u v} }
\end{eqnarray}
のようになりました。総和の前の係数は、弄れるので
\begin{eqnarray}
f \left( v \right) & = & \frac{1}{\sqrt{N}} \sum_{u=0}^{N-1} { F \left( u \right) e^{i \frac{2 \pi}{N} u v} } \\
F \left( u \right) & = & \frac{1}{\sqrt{N}} \sum_{v=0}^{N-1} { f \left( v \right) e^{-i \frac{2 \pi}{N} u v} }
\end{eqnarray}
や
\begin{eqnarray}
f \left( v \right) & = & \frac{1}{N} \sum_{u=0}^{N-1} { F \left( u \right) e^{i \frac{2 \pi}{N} u v} } \\
F \left( u \right) & = & \sum_{v=0}^{N-1} { f \left( v \right) e^{-i \frac{2 \pi}{N} u v} }
\end{eqnarray}
にできます。
計算量の削減
サンプリング数 $N$ が合成数 $\left( N = W \times H \right)$ の場合、$N$ 個の一次元配列を、$W \times H$ の二次元配列と考えます。
F \left( u \right) = \sum_{v=0}^{N-1} { f \left( v \right) e^{-i \frac{2 \pi}{N} u v} }
では、$u$ は $W \times H$ で、 $v$ は $H \times W$ の配列
\begin{eqnarray}
N & = & W H \\
u & = & u_y W + u_x \\
v & = & v_y H + v_x \\
\end{eqnarray}
\\
0 \le \{ u_x, v_y \} \lt W \\
0 \le \{ u_y, v_x \} \lt H
として
\begin{eqnarray}
F \left( u \right) & = & F \left( u_y W + u_x \right) \\
& = & \sum_{v=0}^{N-1} { f \left( v \right) e^{-i \frac{2 \pi}{N} u v} } \\
& = & \sum_{v_x=0}^{H-1} \sum_{v_y=0}^{W-1} { f \left( v_y H + v_x \right) e^{-i \frac{2 \pi}{W H} \left( u_y W + u_x \right) \left( v_y H + v_x \right) } }
\end{eqnarray}
二重総和にします。指数関数の指数部を整理すると
\begin{eqnarray}
& & -i \frac{2 \pi}{W H} \left( u_y W + u_x \right) \left( v_y H + v_x \right) \\
& = & -i \frac{2 \pi}{W H} \left( u_y v_y W H + u_y v_x W + u_x v_y H + u_x v_x \right) \\
& = & -2 \pi i \left( u_y v_y + \frac{u_y v_x}{H} + \frac{u_x v_y}{W} + \frac{u_x v_x}{W H} \right)
\end{eqnarray}
$u_y$ と $v_y$ は整数だから
e^{-2 \pi i \ u_y v_y} = 1
となって（配列が $u$ と $v$ で違うため）省けるので
-2 \pi i \left( \frac{u_y v_x}{H} + \frac{u_x v_y}{W} + \frac{u_x v_x}{W H} \right)
が指数部になります。戻って続けると
\begin{eqnarray}
F \left( u \right) & = & F \left( u_y W + u_x \right) \\
& = & \sum_{v_x=0}^{H-1} \sum_{v_y=0}^{W-1} { f \left( v_y H + v_x \right) e^{-2 \pi i \left( \frac{u_y v_x}{H} + \frac{u_x v_y}{W} + \frac{u_x v_x}{W H} \right)} } \\
& = & \sum_{v_x=0}^{H-1}{ \left[ e^{-2 \pi i \frac{u_x v_x}{W H} } \sum_{v_y=0}^{W-1}{ f \left( v_y H + v_x \right) e^{-2 \pi i \frac{u_x v_y}{W} } } \right] e^{-2 \pi i \frac{u_y v_x}{H} } }
\end{eqnarray}
となりますが、総和の内側を
g \left( u_x, v_x \right) = e^{-2 \pi i \frac{u_x v_x}{W H} } \sum_{v_y=0}^{W-1}{ f \left( v_y H + v_x \right) e^{-2 \pi i \frac{u_x v_y}{W} } }
とすると
F \left( u_y W + u_x \right) = \sum_{v_x=0}^{H-1}{ g \left( u_x, v_x \right) e^{-2 \pi i \frac{u_y v_x}{H} } }
になるので、$g \left( u_x, v_x \right)$ の結果を再利用できます。
サンプリング数 $N$ を素因数分解して $N = P_1 P_2 \cdots P_k$ になる場合、$k$ 重総和にできるため、 $k$ が大きいほど、使い回せる計算結果が多くなって、計算量が減ります。
サンプリング数が 2 の冪乗の場合
サンプリング数が $N = 2^n$ となる場合、$F \left( u \right), f \left( v \right)$ の $u,v$ を
\begin{eqnarray}
u & = & u_0 + 2 u_1 + 2^2 u_2 + \cdots + 2^{n-2} u_{n-2} + 2^{n-1} u_{n-1} \\
v & = & v_0 + 2 v_1 + 2^2 v_2 + \cdots + 2^{n-2} v_{n-2} + 2^{n-1} v_{n-1}
\end{eqnarray}
とします。$u_k,v_k$ は、$0$ か $1$ なので、2進数表現
\begin{eqnarray}
u & = & u_{n-1} u_{n-2} \cdots u_2 u_1 u_0 \\
v & = & v_{n-1} v_{n-2} \cdots v_2 v_1 v_0
\end{eqnarray}
にして簡略化します。離散フーリエ変換は、多重総和
\begin{eqnarray}
F \left( u \right) & = & F \left( u_{n-1} u_{n-2} \cdots u_2 u_1 u_0 \right) \\
& = & \sum_{v_0=0}^1 \sum_{v_1=0}^1 \sum_{v_2=0}^1 \cdots \sum_{v_{n-2}=0}^1 \sum_{v_{n-1}=0}^1 f \left( v_{n-1} v_{n-2} \cdots v_2 v_1 v_0 \right) e^{-2 \pi i \frac{ \left( u_{n-1} u_{n-2} \cdots u_2 u_1 u_0 \right) \left( v_{n-1} v_{n-2} \cdots v_2 v_1 v_0 \right) }{2^n}}
\end{eqnarray}
にすることができます。
少し戻って、二重総和から始めます。$N=W \times H$ への割り当てを
\begin{eqnarray}
n & = & w + h \\
W & = & 2^w \\
H & = & 2^h \\
N & = & W \times H = 2^n = 2^{w+h}
\end{eqnarray}
とします。$u$ は、下位のビット数が $w$ で、上位のビット数が $h$ として、2進数表現と変数を
\begin{eqnarray}
u_H^h & = & u_{n-1} u_{n-2} \cdots u_{w+1} u_w \\
u_L^w & = & u_{w-1} u_{w-2} \cdots u_{2} u_{1} u_{0} \\
u & = & 2^w u_H^h + u_L^w
\end{eqnarray}
と表現することにします。$u$ の添字(上) $w,h$ はビット数で、添字(下) $H,L$ は上位・下位です。$v$ は、上位・下位のビット数を逆にして
\begin{eqnarray}
v_H^w & = & v_{n-1} v_{n-2} \cdots v_{h+1} v_h \\
v_L^h & = & v_{h-1} v_{h-2} \cdots v_{2} v_{1} v_{0} \\
v & = & 2^h v_H^w + v_L^h
\end{eqnarray}
とすると
\begin{eqnarray}
F \left( u \right) & = & F \left( 2^w u_H^h + u_L^w \right) \\
& = & \sum_{v_L^h=0}^{2^h-1}{ \sum_{v_H^w=0}^{2^w-1} { f \left( 2^h v_H^w + v_L^h \right) e^{ -2 \pi i \frac{ \left( 2^w u_H^h + u_L^w \right) \left( 2^h v_H^w + v_L^h \right) }{ 2^n } } } } \\
& = & \sum_{v_L^h=0}^{2^h-1}{ \left[ e^{ -2 \pi i \frac{ u_L^w v_L^h }{ 2^n } } \sum_{v_H^w=0}^{2^w-1} { f \left( 2^h v_H^w + v_L^h \right) e^{ -2 \pi i \frac{ u_L^w v_H^w }{2^w} } } \right] e^{ -2 \pi i \frac{ u_H^h v_L^h }{2^h} } }
\end{eqnarray}
だから、内側の総和を $g_w \left( u_L^w, v_L^h \right)$ として
\begin{eqnarray}
g_w \left( u_L^w, v_L^h \right) & = & e^{ -2 \pi i \frac{ u_L^w v_L^h }{ 2^n } } \sum_{v_H^w=0}^{2^w-1} { f \left( 2^h v_H^w + v_L^h \right) e^{ -2 \pi i \frac{ u_L^w v_H^w }{2^w} } } \\
F \left( u \right) & = & F \left( 2^w u_H^h + u_L^w \right) = \sum_{v_L^h=0}^{2^h-1}{ g_w \left( u_L^w, v_L^h \right) e^{ -2 \pi i \frac{ u_H^h v_L^h }{2^h} } }
\end{eqnarray}
にします。$w=1$ とすると
\begin{eqnarray}
g_1 \left( u_0, v_L^{n-1} \right)
& = & e^{ -2 \pi i \frac{ v_L^{n-1} }{2^n} u_0 } \sum_{ v_{n-1}=0 }^{1}{ f \left( 2^{n-1} v_{n-1} + v_L^{n-1} \right) e^{ -2 \pi i \frac{ v_{n-1} }{2} u_0 } } \\
& = & e^{ -2 \pi i \frac{ v_L^{n-1} }{2^n} u_0 } \left( \ f \left( v_L^{n-1} \right) + f \left( 2^{n-1} + v_L^{n-1} \right) e^{ - \pi i \ u_0 } \right) \\
& = & e^{ -2 \pi i \frac{ v_L^{n-1} }{2^n} u_0 } \left(\ f \left( v_L^{n-1} \right) + \left( -1 \right)^{u_0} f \left( 2^{n-1} + v_L^{n-1} \right) \right)
\end{eqnarray}
となって、$u_0$ に $0$ と $1$ を代入すると
\begin{eqnarray}
g_1 \left( 0, v_L^{n-1} \right) & = & f \left( v_L^{n-1} \right) + f \left( 2^{n-1} + v_L^{n-1} \right) \\
g_1 \left( 1, v_L^{n-1} \right) & = & \left( \ f \left( v_L^{n-1} \right) - f \left( 2^{n-1} + v_L^{n-1} \right) \right) e^{ -\pi i \frac{ v_L^{n-1} }{2^{n-1}} } \\
F \left( u \right) & = & F \left( u_H^{n-1} + u_0 \right) = \sum_{v_L^{n-1}=0}^{2^{n-1}-1}{ g_1 \left( u_0, v_L^{n-1} \right) e^{ -2 \pi i \frac{ u_H^{n-1} v_L^{n-1} }{ 2^n } } }
\end{eqnarray}
を得ます。$w=2$ では
\begin{eqnarray}
u_L^2 & = & u_1 u_0 \\
v_H^2 & = & v_{n-1} v_{n-2} \\
g_2 \left( u_L^2, v_L^{n-2} \right) & = & e^{ -2 \pi i \frac{ u_L^2 v_L^{n-2} }{2^n} } \sum_{ v_H^2=0 }^{3}{ f \left( 2^{n-2} v_H^2 + v_L^{n-2} \right) e^{ -2 \pi i \frac{ u_L^2 v_H^2 }{2^2} } } \\
& = & e^{ -2 \pi i \frac{ u_L^2 v_L^{n-2} }{2^n} } \sum_{ v_{n-2}=0 }^{1}{ \sum_{ v_{n-1}=0 }^{1}{ f \left( 2^{n-2} v_H^2 + v_L^{n-2} \right) e^{ -2 \pi i \frac{ u_L^2 v_H^2 }{2^2} } } }
\end{eqnarray}
ですが、末項の指数関数の指数部の分数の分子は
\begin{eqnarray}
u_L^2 v_H^2
& = & \left( 2 u_1 + u_0 \right) \left( 2 v_{n-1} + v_{n-2} \right) \\
& = & 4 u_1 v_{n-1} + 2 u_0 v_{n-1} + 2 u_1 v_{n-2} + u_0 v_{n-2} \\
& = & 4 u_1 v_{n-1} + 2 u_0 v_{n-1} + u_L^2 v_{n-2}
\end{eqnarray}
これを末項の指数関数に戻すと
\begin{eqnarray}
e^{ -2 \pi i \frac{ u_L^2 v_H^2 }{2^2} }
& = & e^{ -2 \pi i \frac{ 4 u_1 v_{n-1} + 2 u_0 v_{n-1} + u_L^2 v_{n-2} }{2^2} } \\
& = & e^{ -2 \pi i \frac{v_{n-1} }{2} u_0 } e^{ -2 \pi i \frac{ u_L^2 v_{n-2} }{2^2} }
\cdots \left( e^{ -2 \pi i \frac{ 4 u_1 v_{n-1} }{2^2} } = e^{ -2 \pi i \ u_1 v_{n-1} } = 1 \right) \\
g_2 \left( u_L^2, v_L^{n-2} \right) & = & e^{ -2 \pi i \frac{ u_L^2 v_L^{n-2} }{2^n} } \sum_{ v_{n-2}=0 }^{1}{ \left[ e^{ -2 \pi i \frac{ u_L^2 v_{n-2} }{2^2} } \sum_{ v_{n-1}=0 }^{1}{ f \left( 2^{n-2} v_H^2 + v_L^{n-2} \right) e^{ -2 \pi i \frac{ v_{n-1} }{2} u_0 } } \right] }
\end{eqnarray}
となりますが、前方2つの指数関数を
e^{ -2 \pi i \frac{ u_L^2 v_L^{n-2} }{2^n} } e^{ -2 \pi i \frac{ u_L^2 v_{n-2} }{2^2} }
として、指数部を整理すると
\begin{eqnarray}
& & \left( -2 \pi i \frac{ u_L^2 v_L^{n-2} }{2^n} \right) + \left( -2 \pi i \frac{ u_L^2 v_{n-2} }{2^2} \right) \\
& = & -2 \pi i \frac{ 2^{n-2} v_{n-2} + v_L^{n-2} }{2^n} u_L^2 \\
& = & -2 \pi i \frac{ 2^{n-2} v_{n-2} + v_L^{n-2} }{2^n} \left( 2 u_1 + u_0 \right) \\
& = & \left( -2 \pi i \frac{ v_L^{n-2} }{ 2^{n-1} } u_1 \right) + \left( -2 \pi i \frac{ v_{n-2} }{2} u_1 \right) + \left( -2 \pi i \frac{ v_L^{n-1} }{2^n} u_0 \right)
\end{eqnarray}
指数関数に戻して
e^{ -2 \pi i \frac{ u_L^2 v_L^{n-2} }{2^n} } e^{ -2 \pi i \frac{ u_L^2 v_{n-2} }{2^2} } = e^{ -2 \pi i \frac{ v_L^{n-2} }{ 2^{n-1} } u_1 } e^{ -2 \pi i \frac{ v_{n-2} }{2} u_1 } e^{ -2 \pi i \frac{ v_L^{n-1} }{2^n} u_0 }
$g_2 \left( u_L^2, v_L^{n-2} \right)$ に適用すると
\begin{eqnarray}
g_2 \left( u_L^2, v_L^{n-2} \right)
& = & e^{ -2 \pi i \frac{ v_L^{n-2} }{ 2^{n-1} } u_1 } \sum_{ v_{n-2}=0 }^{1}{ \left[ e^{ -2 \pi i \frac{ v_{n-2} }{2} u_1 } \left( e^{ -2 \pi i \frac{ v_L^{n-1} }{2^n} u_0 } \sum_{ v_{n-1}=0 }^{1}{ f \left( 2^{n-2} v_H^2 + v_L^{n-2} \right) e^{ -2 \pi i \frac{ v_{n-1} }{2} u_0 } } \right) \right] } \\
& = & e^{ -2 \pi i \frac{ v_L^{n-2} }{ 2^{n-1} } u_1 } \sum_{ v_{n-2}=0 }^{1}{ g_1 \left( u_0, 2^{n-2} v_{n-2} + v_L^{n-2} \right) e^{ -2 \pi i \frac{ v_{n-2} }{2} u_1 } } \\
& = & e^{ -2 \pi i \frac{ v_L^{n-2} }{ 2^{n-1} } u_1 } \left( g_1 \left( u_0, v_L^{n-2} \right) + g_1 \left( u_0, 2^{n-2} + v_L^{n-2} \right) e^{ - \pi i \ u_1 } \right) \\
& = & e^{ -2 \pi i \frac{ v_L^{n-2} }{ 2^{n-1} } u_1 } \left( g_1 \left( u_0, v_L^{n-2} \right) + \left( -1 \right)^{u_1} g_1 \left( u_0, 2^{n-2} + v_L^{n-2} \right) \right)
\end{eqnarray}
ですが、$u_1$ は、$0$ または $1$ だから
\begin{eqnarray}
g_2 \left( 0 + u_0, v_L^{n-2} \right) & = & g_1 \left( u_0, v_L^{n-2} \right) + g_1 \left( u_0, 2^{n-2} + v_L^{n-2} \right) \\
g_2 \left( 2 + u_0, v_L^{n-2} \right) & = & \left( g_1 \left( u_0, v_L^{n-2} \right) - g_1 \left( u_0, 2^{n-2} + v_L^{n-2} \right) \right) e^{ - \pi i \frac{ v_L^{n-2} }{ 2^{n-2} } } \\
F \left( u \right) = F \left( u_H^{n-2} + u_L^2 \right) & = & \sum_{v_L^{n-2}=0}^{2^{n-2}-1}{ g_2 \left( u_L^2, v_L^{n-2} \right) e^{ -2 \pi i \frac{ u_H^{n-2} v_L^{n-2} }{ 2^n } } }
\end{eqnarray}
を得ます。$g_1$ と $g_2$ を並べると
\begin{eqnarray}
g_1 \left( u_0, v_L^{n-1} \right) & = & e^{ -2 \pi i \frac{ v_L^{n-1} }{2^n} u_0 } \left(\ f \left( v_L^{n-1} \right) + \left( -1 \right)^{u_0} f \left( 2^{n-1} + v_L^{n-1} \right) \right) \\
g_2 \left( u_L^2, v_L^{n-2} \right) & = & e^{ -2 \pi i \frac{ v_L^{n-2} }{ 2^{n-1} } u_1 } \left( g_1 \left( u_0, v_L^{n-2} \right) + \left( -1 \right)^{u_1} g_1 \left( u_0, 2^{n-2} + v_L^{n-2} \right) \right)
\end{eqnarray}
ここで、$g_0 \left( 0, v \right) = f \left( v \right)$ とすると
\begin{eqnarray}
g_1 \left( u_0, v_L^{n-1} \right) & = & e^{ -2 \pi i \frac{ v_L^{n-1} }{2^n} u_0 } \left( g_0 \left( 0, v_L^{n-1} \right) + \left( -1 \right)^{u_0} g_0 \left( 0, 2^{n-1} + v_L^{n-1} \right) \right) \\
g_2 \left( u_L^2, v_L^{n-2} \right) & = & e^{ -2 \pi i \frac{ v_L^{n-2} }{ 2^{n-1} } u_1 } \left( g_1 \left( u_0, v_L^{n-2} \right) + \left( -1 \right)^{u_1} g_1 \left( u_0, 2^{n-2} + v_L^{n-2} \right) \right)
\end{eqnarray}
とできるので、$g_s$ は
g_s \left( u_L^s, v_L^{n-s} \right) = e^{ -2 \pi i \frac{ v_L^{n-s} }{ 2^{n - \left( s-1\right)} } u_{s-1} } \left( g_{s-1} \left( u_L^{s-1}, v_L^{n - \left( s-1 \right)} \right) + \left( -1 \right)^{u_{s-1}} g_{s-1} \left( u_L^{s-1}, 2^{n-s} + v_L^{n - \left( s-1 \right)} \right) \right)
になります。
$u,v$ を2進数として、サンプル数が16個の場合は
\begin{eqnarray}
u & = & u_3 u_2 u_1 u_0 \\
v & = & v_3 v_2 v_1 v_0 \\
\\
g_1 \left( 0, v_2 v_1 v_0 \right) & = & f \left( 0 v_2 v_1 v_0 \right) + f \left( 1 v_2 v_1 v_0 \right) \\
g_1 \left( 1, v_2 v_1 v_0 \right) & = & \left( g_0 \left( 0, 0 v_2 v_1 v_0 \right) - g_0 \left( 0, 1 v_2 v_1 v_0 \right) \right) e^{ - \pi i \frac{ v_2 v_1 v_0 }{2^3} } \\
\\
g_2 \left( 0 u_0, v_1 v_0 \right) & = & \left( g_1 \left( u_0, 0 v_1 v_0 \right) + g_1 \left( u_0, 1 v_1 v_0 \right) \right) \\
g_2 \left( 1 u_0, v_1 v_0 \right) & = & \left( g_1 \left( u_0, 0 v_1 v_0 \right) - g_1 \left( u_0, 1 v_1 v_0 \right) \right) e^{ - \pi i \frac{ v_1 v_0 }{2^2} } \\
\\
g_3 \left( 0 u_1 u_0, v_0 \right) & = & \left( g_2 \left( u_1 u_0, 0 v_0 \right) + g_2 \left( u_1 u_0, 1 v_0 \right) \right) \\
g_3 \left( 1 u_1 u_0, v_0 \right) & = & \left( g_2 \left( u_1 u_0, 0 v_0 \right) - g_2 \left( u_1 u_0, 1 v_0 \right) \right) e^{ - \pi i \frac{ v_0 }{2^1} } \\
\\
g_4 \left( 0 u_2 u_1 u_0 \right) & = & g_3 \left( u_2 u_1 u_0, 0 \right) + g_3 \left( u_2 u_1 u_0, 1 \right) \\
g_4 \left( 1 u_2 u_1 u_0 \right) & = & g_3 \left( u_2 u_1 u_0, 0 \right) - g_3 \left( u_2 u_1 u_0, 1 \right) \\
\\
F \left(u_3 u_2 u_1 u_0 \right) & = & g_4 \left(u_3 u_2 u_1 u_0 \right)
\end{eqnarray}
プログラム風に表現すると
# uuuu や vvvv などは、2進数表現
for vvv in range(8):
g1(0,vvv) = f(0vvv) + f(1vvv)
g1(1,vvv) = (f(0vvv) - f(1vvv)) * Exp(- PI * I * vvv / 8)
for uvv in range(8):
g2(0u,vv) = g1(u,0vv) + g1(u,1vv)
g2(1u,vv) = (g1(u,0vv) - g1(u,1vv)) * Exp(- PI * I * vv / 4)
for uuv in range(8):
g3(0uu,v) = g2(uu,0v) + g2(uu,1v)
g3(1uu,v) = (g2(uu,0v) - g2(uu,1v)) * Exp(- PI * I * v / 2)
for uuu in range(8):
g4(0uuu) = g3(uuu,0) + g3(uuu,1)
g4(1uuu) = (g3(uuu,0) - g3(uuu,1)) # * (Exp(-PI * I * 0 / 1) = 1)
for uuuu in range(16):
F(uuuu) = g4(uuuu)
になり、$f$ から $g_1$ を(8回ループで)、$g_1$ から $g_2$ を(8回ループで)、$\cdots$ という順番で計算します。
サンプル プログラム
サンプル プログラムを Python3 で記述してみる。
test.py
#!/usr/bin/env python3
import cmath
def to_compliex_list(data):
"""様々な内容(型)のリストを複素数表にする"""
r = []
for p in data:
if type(p) == complex:
r.append(p)
elif type(p) in (int, float):
r.append(complex(float(p)))
elif len(p) < 2:
r.append(complex(float(p[0])))
else:
r.append(complex(float(p[0]), float(p[1])))
return r
def dft(data, inv=False):
"""離散フーリエ変換/逆変換"""
N = len(data)
M = 1.0 if inv else 1 / N
s2_pi = 2 * cmath.pi * (1.0 if inv else -1.0)
if False:
return [M * sum(d * cmath.exp(complex(0, s2_pi * ((u * v) % N) / N))
for v, d in enumerate(data)) for u in range(N)]
r = []
for u in range(N):
f = complex()
for v in range(N):
p = (u * v) % N
e = complex(0, s2_pi * p / N)
f += data[v] * cmath.exp(e)
r.append(f * M)
return r
def fft_p2(data, inv=False):
"""高速離散フーリエ変換/逆変換"""
N = len(data)
if (N & -N) != N:
raise ValueError(f'データ長({N})が 2 の冪乗ではありません')
NB = N.bit_length() - 1
s_pi = cmath.pi * (1.0 if inv else -1.0)
coeff = 1.0 if inv else 1.0 / N
dp = [d * coeff for d in data]
NS = N >> 1
NU = 1
NV = N >> 1
for g in range(NB):
dn = [None] * N
for u in range(NU):
un = u * NV
up = un << 1
for v in range(NV):
vp = up + v
vn = un + v
p0 = dp[vp]
p1 = dp[vp + NV]
n0 = (p0 + p1)
n1 = (p0 - p1) * cmath.exp(complex(0, s_pi * v / NV))
dn[vn] = n0
dn[vn + NS] = n1
dp = dn
NU <<= 1
NV >>= 1
return dn
def dump(data, msg=''):
if len(msg):
print(msg)
for n in range(len(data)):
print('%3d: %+10.8f, %+10.8f' % (n, data[n].real, data[n].imag))
def test():
N = 16
l = [- 1.0 + 2 * float(n) / N for n in range(N)]
dat = to_compliex_list(l)
dump(dat, 'DATA')
dfn = dft(dat, False)
dump(dfn, 'DFT')
ffn = fft_p2(dat, False)
dump(ffn, 'FFT')
dfi = dft(dfn, True)
dump(dfi, 'IDFT')
ffi = fft_p2(ffn, True)
dump(ffi, 'IFFT')
if __name__ == '__main__':
test()
実行結果
$ python3 test.py
DATA
0: -1.00000000, +0.00000000
1: -0.87500000, +0.00000000
2: -0.75000000, +0.00000000
3: -0.62500000, +0.00000000
4: -0.50000000, +0.00000000
5: -0.37500000, +0.00000000
6: -0.25000000, +0.00000000
7: -0.12500000, +0.00000000
8: +0.00000000, +0.00000000
9: +0.12500000, +0.00000000
10: +0.25000000, +0.00000000
11: +0.37500000, +0.00000000
12: +0.50000000, +0.00000000
13: +0.62500000, +0.00000000
14: +0.75000000, +0.00000000
15: +0.87500000, +0.00000000
DFT
0: -0.06250000, +0.00000000
1: -0.06250000, +0.31420872
2: -0.06250000, +0.15088835
3: -0.06250000, +0.09353786
4: -0.06250000, +0.06250000
5: -0.06250000, +0.04176116
6: -0.06250000, +0.02588835
7: -0.06250000, +0.01243202
8: -0.06250000, -0.00000000
9: -0.06250000, -0.01243202
10: -0.06250000, -0.02588835
11: -0.06250000, -0.04176116
12: -0.06250000, -0.06250000
13: -0.06250000, -0.09353786
14: -0.06250000, -0.15088835
15: -0.06250000, -0.31420872
FFT
0: -0.06250000, +0.00000000
1: -0.06250000, +0.31420872
2: -0.06250000, +0.15088835
3: -0.06250000, +0.09353786
4: -0.06250000, +0.06250000
5: -0.06250000, +0.04176116
6: -0.06250000, +0.02588835
7: -0.06250000, +0.01243202
8: -0.06250000, +0.00000000
9: -0.06250000, -0.01243202
10: -0.06250000, -0.02588835
11: -0.06250000, -0.04176116
12: -0.06250000, -0.06250000
13: -0.06250000, -0.09353786
14: -0.06250000, -0.15088835
15: -0.06250000, -0.31420872
IDFT
0: -1.00000000, -0.00000000
1: -0.87500000, +0.00000000
2: -0.75000000, -0.00000000
3: -0.62500000, -0.00000000
4: -0.50000000, -0.00000000
5: -0.37500000, +0.00000000
6: -0.25000000, +0.00000000
7: -0.12500000, +0.00000000
8: -0.00000000, +0.00000000
9: +0.12500000, +0.00000000
10: +0.25000000, +0.00000000
11: +0.37500000, -0.00000000
12: +0.50000000, +0.00000000
13: +0.62500000, +0.00000000
14: +0.75000000, -0.00000000
15: +0.87500000, -0.00000000
IFFT
0: -1.00000000, +0.00000000
1: -0.87500000, +0.00000000
2: -0.75000000, +0.00000000
3: -0.62500000, +0.00000000
4: -0.50000000, -0.00000000
5: -0.37500000, +0.00000000
6: -0.25000000, +0.00000000
7: -0.12500000, +0.00000000
8: -0.00000000, +0.00000000
9: +0.12500000, -0.00000000
10: +0.25000000, -0.00000000
11: +0.37500000, +0.00000000
12: +0.50000000, +0.00000000
13: +0.62500000, -0.00000000
14: +0.75000000, -0.00000000
15: +0.87500000, -0.00000000
DFT と FFT が同じになったので、よさそうです。