円周率.jp (http://xn--w6q13e505b.jp/method/fft/complexity.html)

計算量の導出

ここでは FFT の計算量が $O(n\log n)$ になることについての簡単な説明を行う.

ちなみに, 多倍長乗算はフーリエ変換と乗算との関係に記載したように,

  1. 2 回の FFT $\cdots O(n\log n)$
  2. 要素ごとの積 $\cdots O(n)$
  3. 逆 FFT $\cdots O(n\log n)$

という手順で行うことができるので,その計算量は $O(n\log n)$ になる.

この計算量を示す過程には 2 種類の導き方があるので, その両者について説明する. その際,変換元のデータは時間軸上のデータ, 変換後のデータは周波数軸上のデータ, というように捉えると名前が分かりやすい.

時間間引きFFT

定義式のうち,$j$ が奇数のものと偶数のものに分ける.

\[ A_k = \sum_{j=0}^{n/2-1} a_{2j}\omega^{2jk} + \omega^k\sum_{j=0}^{n/2-1}a_{2j+1}\omega^{2jk} \tag{1} \]

こうなると右辺の $\sum$ の中は項数 $n/2$ の DFT になっているので,$A_k$ の方もそれに合わせる. ただしこちら側も $k$ の偶奇で分割しては意味がないので $k$ と $k+n/2$ に分けて $0\leq k \lt n/2$ にする.

\[ A_{k+n/2} = \sum_{j=0}^{n/2-1} a_{2j}\omega^{2jk} - \omega^k\sum_{j=0}^{n/2-1}a_{2j+1}\omega^{2jk} \tag{2} \]

ここで (1) 式と (2) 式を同時に求めることを考えると, ${\displaystyle \sum_{j=0}^{n/2-1}} a_{2j}\omega^{2jk}$ と $\omega^k {\displaystyle \sum_{j=0}^{n/2-1}} a_{2j+1}\omega^{2jk}$ を求めればよいことが分かる. また,それらは項数 $n/2$ の DFT であるため, この方法を再帰的に適用することができることも分かる. プログラムで書くと

// x[n] の DFT を計算する時間間引き FFT
void InTimeFFT(int n, Complex* x, Complex* y) {
  if (n <= 1) return;

  Complex *x0 = x, *x1 = &x[n/2];
  Complex *y0 = y, *y1 = &y[n/2];
  for (int j = 0; j < n / 2; ++j) {
    y0[j] = x[2 * j];
    y1[j] = x[2 * j + 1];
  }
  InTimeFFT(n / 2, y0, x0);
  InTimeFFT(n / 2, y1, x1);
  for (int k = 0; k < n / 2; ++k) {
    double wt = -2 * M_PI * k / n;
    Complex t = y1[k] * Complex(cos(wt), sin(wt));
    x0[k] = y0[k] + t;
    x1[k] = y0[k] - t;
  }
}

である。y[n] は一時領域として使う。

ここで項数 $n$ の FFT の計算量を $T(n)$ とすると,先ほどの再帰的な FFT 分割適用は

\[ T(n) = 2T(n/2) + Cn \]

となる.末尾の $Cn$ は $\omega^k$ の乗算や 2 つの FFT の結果の加減算にかかるコストを示している. この式を解くと

\[ T(n) = O(n\log n) \]

となる.

周波数間引きFFT

今度は逆に,定義式のうち $j$ を $j$ と $j+n/2$ とに分け,$k$ を偶奇に分けてみる.

\[ \begin{eqnarray} A_{2k} &=& {\displaystyle \sum_{j=0}^{n/2-1}} a_j\omega^{2jk} &+& {\displaystyle \sum_{j=0}^{n/2-1}} a_{j+n/2}\omega^{2jk} &=& {\displaystyle \sum_{j=0}^{n/2-1}} \left( a_j + a_{j+n/2} \right) \omega^{2jk} \tag{3}\\ A_{2k+1} &=& {\displaystyle \sum_{j=0}^{n/2-1}} a_j\omega^j\omega^{2jk} &-& {\displaystyle \sum_{j=0}^{n/2-1}} a_{j+n/2}\omega^j\omega^{2jk} &=& {\displaystyle \sum_{j=0}^{n/2-1}} \left( a_j - a_{j+n/2} \right) \omega^j \omega^{2jk} \tag{4} \end{eqnarray} \]

ここで (3) 式と (4) 式については

\[ \begin{cases} a_j^{(0)} = a_j + a_{j+n/2}\\ a_j^{(1)} = a_j - a_{j+n/2} \omega^j \end{cases} \]

の演算は同時にできる. これをプログラムコードにすると

// x[n] の DFT を計算する周波数間引き FFT
void InFreqFFT(int n, Complex* x, Complex* y) {
  if (n <= 1)
    return;

  Complex *y0 = y, *y1 = &y[n/2];
  for (int j = 0; j < n / 2; ++j) {
    double wt = -2 * M_PI * j / n;
    y0[j] = x0[j] + x1[j];
    y1[j] = (x0[j] - x1[j]) * Complex(cos(wt), sin(wt));
  }
  InFreqFFT(n / 2, y0, x0);
  InFreqFFT(n / 2, y1, x1);
  for (int k = 0; k < n / 2; ++k) {
    x[2*k  ] = y0[k];
    x[2*k+1] = y1[k];
  }
}

となる.