メモの日々


2024年04月01日(月) [長年日記]

[c++][dev][math] 高速フーリエ変換

高速フーリエ変換についてメモ。

高速フーリエ変換とは、離散フーリエ変換を高速に計算するアルゴリズムのこと。

離散フーリエ変換

離散フーリエ変換とは、項数 N の複素数列 \((a_n)_{n=0}^{N-1}\) から、次で定める項数 N の複素数列 \((b_n)_{n=0}^{N-1}\) を求める操作のこととする。

\[ \ b_n = \sum_{x=0}^{N-1}{a_x \exp(-i \frac{2 \pi nx}{N})} \]

\(W_N=\exp\left(-i\frac{2 \pi}{N}\right)\) とおくと、この式は次のように書ける。

\[ \ b_n = \sum_{x=0}^{N-1}{W_N^{nx} a_x} \]

この式の掛け算の回数は N 回で、すべての 0≦n≦N に対して計算するには N×N 回の掛け算が必要になる。

高速フーリエ変換の解説は N を2のべき乗に制限するものが多いが、ここでは非素数の N に対して(素朴に計算するのに比べて)高速に計算できるアルゴリズムを考える。なお、以下の話は N が素数であっても適用できる(R=N、S=1 とすればよい)が、その場合は計算量が減らない。

計算量を減らすための式変形

R を N の約数、S を N/R とする。求めたい複素数列 \((b_n)_{n=0}^{N-1}\) を項数 S の R 個の複素数列 \((b_{Rm})_{m=0}^{S-1}, (b_{Rm+1})_{m=0}^{S-1}, ..., (b_{Rm+R-1})_{m=0}^{S-1} \) に分割する。数列 \((b_{Rm+k})_{m=0}^{S-1}\) の一般項は次のように表せる(この式変形の導出は本稿の最後に示す)。

\[ \ b_{Rm+k} = \sum_{y=0}^{S-1}{W_S^{my} \left( W_N^{ky} \sum_{x=0}^{R-1}{W_R^{kx} a_{Sx+y}} \right)} \]

この式の掛け算の回数を考えると、

  • 括弧内の計算に R+1 回の掛け算が必要で、これをすべての 0≦y≦S-1 と 0≦k≦R-1 に対して計算するので全体では (R+1)×S×R 回
  • \(W_S^{my}\)との積が S 回でこれをすべての 0≦m≦S-1 と 0≦k≦R-1 に対して計算するので全体では S×S×R 回

になる(と思う)。合わせて (R+1)×S×R + S×S×R = (R+S+1)×N 回の掛け算が必要で、これはほとんどのケースで元の式の N×N 回より少ない回数になる。括弧内の計算がmに依存しておらず計算結果を再利用できることが効いている。ただし、元の式には無かった \(W_S^{my}\) と \(W_R^{kx}\) の計算が新たに発生していることには注意。また、\((b_n)_{n=0}^{N-1}\) を分割する際に並び順が変わっているので、並び替えのコストも発生する。

ここで、\(b_{Rm+k}\) の式が項数 S の複素数列 \(\left( W_N^{ky} \sum_{x=0}^{R-1}{W_R^{kx} a_{Sx+y}} \right)_{y=0}^{S-1}\) に対する離散フーリエ変換の形になっていることに注目してほしい。これにより、R個の数列 \((b_{Rm+k})_{m=0}^{S-1}\) それぞれの計算に対しても上記と同じ式変形で計算量を減らすことができ、この操作は数列の要素数が1になるまで再帰的に行えることになる。

高速フーリエ変換の実装

上のようにして計算量を減らしたコードと減らさず素朴に計算するコードの比較を行うプログラムを書いた。200行くらいあるので全体はGistに置きここにコードを抜粋する。コンパイルにはC++20が必要だが、もっと古いC++でコンパイルできるように直すのは難しくない。

高速フーリエ変換の計算をする関数 fft() は次の通り。

// inputの値を離散フーリエ変換した結果を返す
template<typename T>
std::vector<std::complex<T>> fft(
    std::vector<std::complex<T>> input,
    const twiddle_factor_map<T>& tf_map)
{
    std::vector<std::complex<double>> work(input.size());
    detail::fft(
        input.begin(),
        input.end(),
        work.begin(),
        detail::find_divisor(input.size()),
        tf_map);
    return input;
}
  • 入力をvectorで与え結果をvectorで受け取る。
  • \( W_N \) のべき乗については事前に計算したものを tf_map というマップに設定して渡すようにした。
  • 計算の作業領域用に入力と同じサイズのvectorを使用している。

上で呼び出している実際の処理を行う関数 detail::fft() は次の通り。

// 範囲[first, last)の値を離散フーリエ変換した結果を範囲[first, last)に上書きする。
// work_firstは[first, last)と同じサイズの作業用範囲の先頭。
// radixには範囲[first, last)の要素数の約数を指定する。
template<typename T, typename Iter>
void fft(
    Iter first,
    Iter last,
    Iter work_first,
    std::size_t radix,
    const twiddle_factor_map<T>& tf_map)
{
    const auto n = static_cast<std::size_t>(std::distance(first, last));
    if (n <= 1) return;

    const auto& tfs_r = tf_map.at(radix);
    const auto& tfs_n = tf_map.at(n);
    const auto s = n / radix;

    for (std::size_t i = 0; i < s; ++i) {
        auto it = work_first + i;
        for (std::size_t j = 0; j < radix; ++j) {
            *it = 0;
            for (std::size_t k = 0; k < radix; ++k) {
                *it += *(first + i + s * k) * tfs_r[j * k % radix];
            }
            *it *= tfs_n[i * j];
            it += s;
        }
    }

    for (std::size_t i = 0; i < radix; ++i) {
        const auto it = work_first + s * i;
        fft(it, it + s, first, find_divisor(s), tf_map);
    }

    for (std::size_t i = 0; i < s; ++i) {
        auto it = first + radix * i;
        for (std::size_t j = 0; j < radix; ++j) {
            *it = *(work_first + i + s * j);
            ++it;
        }
    }
}
  • detail::fft() は自分自身を呼び出している。再帰呼び出しの中でvectorを作り直さなくていいように、イテレータ範囲でデータを渡すようにしている。
  • 作業領域内にR個の複素数列 \(\left( W_N^{ky} \sum_{x=0}^{R-1}{W_R^{kx} a_{Sx+y}} \right)_{y=0}^{S-1}\) を設定している。
  • 再帰呼び出しでは作業領域を入力データ、入力データ範囲を作業領域として渡している。
  • 再帰呼び出しが終わったら作業領域内のデータを並び替えたものを計算結果として入力データ範囲に設定している。

実行結果

Gistに置いたファイルの実行結果は以下のような感じ。素朴に計算した場合の実行時間(dft)を高速に実行した場合の実行時間(fft)で割った値を ratio として出力している。delta は計算結果の差の絶対値の最大値。

なお、NumPyのFFTルーチンと比べると遅い(項数 129140163 = 3の17乗 のときに上のコードは23秒、NumPyは10秒で計算)ので、上のコードもしくは本アルゴリズムは改善の余地が大いにある。

入力の項数が少ないうちは素朴に計算する方が速い。

size    2       ratio      0.379        dft        0.000        fft        0.000        delta   0.000e+00
size    3       ratio      0.858        dft        0.000        fft        0.000        delta   0.000e+00
size    4       ratio      0.897        dft        0.000        fft        0.000        delta   0.000e+00
size    5       ratio      0.950        dft        0.000        fft        0.000        delta   0.000e+00
size    6       ratio      0.400        dft        0.000        fft        0.000        delta   7.589e-15
size    7       ratio      0.907        dft        0.000        fft        0.000        delta   0.000e+00
size    8       ratio      0.640        dft        0.000        fft        0.000        delta   3.794e-15
size    9       ratio      2.072        dft        0.000        fft        0.000        delta   0.000e+00
size    10      ratio      0.621        dft        0.000        fft        0.000        delta   7.944e-15
size    11      ratio      0.966        dft        0.000        fft        0.000        delta   0.000e+00
size    12      ratio      0.600        dft        0.000        fft        0.000        delta   2.197e-14
size    13      ratio      1.021        dft        0.000        fft        0.000        delta   0.000e+00
size    14      ratio      0.704        dft        0.000        fft        0.000        delta   1.507e-14
size    15      ratio      0.524        dft        0.000        fft        0.000        delta   1.888e-14

項数2000くらいで、大きくて40倍くらいの差が出てくる。項数が素数の場合は速くならない。

size    1987    ratio      0.998        dft        0.010        fft        0.010        delta   0.000e+00
size    1988    ratio     22.035        dft        0.010        fft        0.000        delta   3.122e-12
size    1989    ratio     46.875        dft        0.010        fft        0.000        delta   2.890e-12
size    1990    ratio      8.610        dft        0.010        fft        0.001        delta   2.621e-12
size    1991    ratio     10.168        dft        0.009        fft        0.001        delta   2.390e-12
size    1992    ratio     19.410        dft        0.010        fft        0.000        delta   2.912e-12
size    1993    ratio      1.000        dft        0.009        fft        0.009        delta   0.000e+00
size    1994    ratio      1.997        dft        0.009        fft        0.005        delta   3.429e-12
size    1995    ratio     46.528        dft        0.010        fft        0.000        delta   2.430e-12
size    1996    ratio      3.932        dft        0.010        fft        0.002        delta   4.590e-12
size    1997    ratio      1.039        dft        0.010        fft        0.009        delta   0.000e+00
size    1998    ratio     33.008        dft        0.010        fft        0.000        delta   2.993e-12
size    1999    ratio      0.892        dft        0.010        fft        0.011        delta   0.000e+00

項数が2のべき乗の場合。数が大きくなるとどんどん差がつく。

--- 2のべき乗 ---
size    2       ratio      1.911        dft        0.000        fft        0.000        delta   0.000e+00
size    4       ratio      0.982        dft        0.000        fft        0.000        delta   0.000e+00
size    8       ratio      0.741        dft        0.000        fft        0.000        delta   7.105e-15
size    16      ratio      0.908        dft        0.000        fft        0.000        delta   1.137e-14
size    32      ratio      1.188        dft        0.000        fft        0.000        delta   6.071e-14
size    64      ratio      2.473        dft        0.000        fft        0.000        delta   2.050e-13
size    128     ratio      4.226        dft        0.000        fft        0.000        delta   3.422e-13
size    256     ratio      8.059        dft        0.000        fft        0.000        delta   1.833e-12
size    512     ratio     15.466        dft        0.001        fft        0.000        delta   2.854e-11
size    1024    ratio     28.213        dft        0.003        fft        0.000        delta   8.084e-11
size    2048    ratio     56.931        dft        0.012        fft        0.000        delta   3.748e-10
size    4096    ratio    104.531        dft        0.044        fft        0.000        delta   7.867e-10
size    8192    ratio    204.091        dft        0.172        fft        0.001        delta   1.531e-09
size    16384   ratio    361.594        dft        0.644        fft        0.002        delta   2.899e-08
size    32768   ratio    715.049        dft        2.648        fft        0.004        delta   6.994e-08
size    65536   ratio   1349.311        dft       11.721        fft        0.009        delta   3.444e-07

項数が3のべき乗の場合。

--- 3のべき乗 ---
size    3       ratio      1.306        dft        0.000        fft        0.000        delta   0.000e+00
size    9       ratio      1.336        dft        0.000        fft        0.000        delta   0.000e+00
size    27      ratio      1.078        dft        0.000        fft        0.000        delta   4.550e-14
size    81      ratio      2.865        dft        0.000        fft        0.000        delta   2.558e-13
size    243     ratio      8.203        dft        0.000        fft        0.000        delta   2.752e-12
size    729     ratio     23.335        dft        0.001        fft        0.000        delta   1.347e-11
size    2187    ratio     57.941        dft        0.013        fft        0.000        delta   1.113e-10
size    6561    ratio    169.265        dft        0.116        fft        0.001        delta   2.806e-09
size    19683   ratio    466.513        dft        1.041        fft        0.002        delta   4.331e-08
size    59049   ratio   1307.727        dft        9.614        fft        0.007        delta   8.864e-08

項数の素因数が異なる数になる場合。

--- 素数の積 ---
size    2       ratio      1.211        dft        0.000        fft        0.000        delta   0.000e+00
size    6       ratio      0.452        dft        0.000        fft        0.000        delta   7.119e-15
size    30      ratio      1.100        dft        0.000        fft        0.000        delta   3.197e-14
size    210     ratio      6.740        dft        0.000        fft        0.000        delta   2.975e-12
size    2310    ratio     61.381        dft        0.014        fft        0.000        delta   6.370e-11
size    30030   ratio    576.535        dft        2.444        fft        0.004        delta   1.515e-07

項数が素数の場合は同じくらいの速さ。

--- 素数 ---
size    2       ratio      0.943        dft        0.000        fft        0.000        delta   0.000e+00
size    3       ratio      0.777        dft        0.000        fft        0.000        delta   0.000e+00
size    5       ratio      0.982        dft        0.000        fft        0.000        delta   0.000e+00
size    7       ratio      1.054        dft        0.000        fft        0.000        delta   0.000e+00
size    11      ratio      1.069        dft        0.000        fft        0.000        delta   0.000e+00
size    13      ratio      1.436        dft        0.000        fft        0.000        delta   0.000e+00
size    1009    ratio      0.999        dft        0.003        fft        0.003        delta   0.000e+00
size    10009   ratio      0.967        dft        0.269        fft        0.278        delta   0.000e+00
size    100003  ratio      0.988        dft       28.695        fft       29.054        delta   0.000e+00

参考

次のページが参考になった。

上で行った変換は混合基数アルゴリズムに該当する。

式変形の詳細

最後に、先に示した \((b_{Rm+k})_{m=0}^{S-1}\) の一般項の導出をメモしておく。

計算の過程で次を使う。xは整数、Nは自然数、RはNの約数、Sは N/R とする。

\[ \ W_N^{Nx} = 1, \quad \ W_N^{Rx} = W_S^x, \quad \ \sum_{x=0}^{N-1}{f(x)} = \sum_{y=0}^{S-1}{\sum_{x=0}^{R-1}{f(Sx+y)}} \]

これらを使うことで次のように変形できる。

\begin{align} \ b_{Rm+k} &= \sum_{x=0}^{N-1}{W_N^{(Rm+k)x} a_x} \\ \ &= \sum_{y=0}^{S-1}{\sum_{x=0}^{R-1}{W_N^{(Rm+k)(Sx+y)} a_{Sx+y}}} \\ \ &= \sum_{y=0}^{S-1}{\sum_{x=0}^{R-1}{W_N^{Nmx} W_N^{Rmy} W_N^{kSx} W_N^{ky} a_{Sx+y}}} \\ \ &= \sum_{y=0}^{S-1}{\sum_{x=0}^{R-1}{W_S^{my} W_R^{kx} W_N^{ky} a_{Sx+y}}} \\ \ &= \sum_{y=0}^{S-1}{W_S^{my} W_N^{ky} \sum_{x=0}^{R-1}{W_R^{kx} a_{Sx+y}}} \\ \end{align}