備忘録 : 高速フーリエ変換

「FFT ってどんなんだったっけ？」となったときに思い出すための備忘録。

# 偶関数と奇関数

偶関数は

```math
f(x) = f(-x)
```

となる関数。

```math
\cos{x} = \cos{\left(-x\right)}
```

奇関数は

```math
\displaylines {
-f(x) = f(-x) \\
\int_{-L}^L f \left( \theta \right) \ d \theta = 0
}
```

となる関数。

```math
- \sin{x} = \sin{\left(-x\right)}
```

偶関数と奇関数の積は奇関数になるので

```math
\displaylines {
\int_{-\pi}^{\pi} \cos{m \theta} \sin{n \theta} \ d \theta = 0 \\
\left( m,n は整数とする \right)
}
```

偶関数と偶関数の積は偶関数で

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

奇関数と奇関数の積は偶関数で

```math
\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)$ の範囲で

```math
f \left( \theta \right) = C + \sum_{n=1}^{\infty} \left( a_n \cos{n \theta} + b_n \sin{n \theta} \right)
```

にできる周期関数として

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

を求めると

```math
\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))$以外はゼロになります。同様に

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

では

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

ゼロの項を省いて整理すると

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

だから

```math
\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$ を負の方向に拡張して

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

これを後で使います。

三角関数を指数関数表現

```math
\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)$ に代入すると

```math
\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)$ を使うと

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

とできます。さらに

```math
c_n = \frac{a_n - i\ b_n}{2}
```

とすると

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

だから

```math
\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)$ にするため

```math
\begin{eqnarray}
x & = & \frac{L}{\pi} \theta \\
\theta & = & \frac{\pi}{L} x \\
d \theta & = & \frac{\pi}{L} dx
\end{eqnarray}
```

により

```math
\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)$ を

```math
f_L \left( x \right) → f \left( x \right) \cdots \left( -L \le x \le L \right)
```

に置き換えて

```math
\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)$ にして、周期性のない関数も対象にします。

```math
\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$ } に変更して

```math
\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$ とすると

```math
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$ では連続したグラフになります。

```math
\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$ を展開すると

```math
\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$ として整理すると

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

内側の積分がフーリエ変換

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

外側の積分が逆フーリエ変換

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

になります。内外の違いは指数関数の指数部の符号だけ。

# 離散フーリエ変換

サンプリング データをフーリエ変換する。

## 基本計算

フーリエ級数まで戻って

```math
\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$ とすると

```math
\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)$ とすると

```math
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)$ は

```math
f \left( x \right) = \sum_{n=0}^{N-1} { c_n e^{i \frac{2 \pi}{N} n x} }
```

となります。これでフーリエ変換と逆変換

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

のようになりました。総和の前の係数を弄って

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

や

```math
\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$ の二次元配列と考えます。

```math
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$ の配列

```math
\begin{eqnarray}
N & = & W H \\
u & = & u_y W + u_x \\
v & = & v_y H + v_x \\
\end{eqnarray}
```

```math:
\begin{eqnarray}
0 \le u_x \lt W \ & , & \ 0 & \le & u_y \lt H \\
0 \le v_x \lt H \ & , & \ 0 & \le & v_y \lt W \\
\end{eqnarray}
```

として

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

二重総和にします。指数関数の指数部を整理すると

```math
\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$ は整数だから

```math
e^{-2 \pi i \ u_y v_y} = 1
```

となって（配列が $u$ と $v$ で違うため）省けるので

```math
-2 \pi i \left( \frac{u_y v_x}{H} + \frac{u_x v_y}{W} + \frac{u_x v_x}{W H} \right)
```

が指数部になります。戻って続けると

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

となりますが、総和の内側を

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

とすると

```math
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$ を

```math
\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進数表現

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

にして簡略化します。離散フーリエ変換は多重総和

```math
\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$ への割り当てを

```math
\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進数表現と変数を

```math
\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$ は、上位・下位のビット数を逆にして

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

とすると

```math
\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)$ として

```math
\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$ とすると

```math
\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$ を代入すると

```math
\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$ では

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

ですが、末項の指数関数の指数部の分数の分子は

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

これを末項の指数関数に戻すと

```math
\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つの指数関数を

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

として、指数部を整理すると

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

指数関数に戻して

```math
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)$ に適用すると

```math
\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$ だから

```math
\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$ を並べると

```math
\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)$ とすると

```math
\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$ は

```math
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個の場合は

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

プログラム風に表現すると

```py

# 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 で記述してみる。

```py: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 が同じになったので、よさそうです。