高速フーリエ変換

高速フーリエ変換#

高速フーリエ変換は離散フーリエ変換の高速なアルゴリズム

1の原始 $N$ 乗根による表現#

1の原始 $N$ 乗根

$ω_{N}$ を次のように定義する。

ω_{N} := e^{i 2 π / N}

これは $N$ 次方程式 $z^{N} - 1 = 0$ の解であり、次の関係を満たす。

ω_{N}^{N} = 1, ω_{N}^{k} \neq 1, k = 1, 2, \dots, N - 1

上の式を満たす $ω_{N}$ を 1の原始 $N$ 乗根 と呼ぶ。

$ω_{N}$ を用いると、離散フーリエ変換の式

f_{l} = \sum_{k = 0}^{N - 1} F_{k} e^{i 2 π k l / N}, F_{k} = \frac{1}{N} \sum_{l = 0}^{N - 1} f_{l} e^{- i 2 π k l / N}

は次のように書くことができる

f_{l} = \sum_{k = 0}^{N - 1} F_{k} ω_{N}^{k l}, F_{k} = \frac{1}{N} \sum_{l = 0}^{N - 1} f_{l} ω_{N}^{- k l}

Tip

フーリエ変換は線形変換

上記だけだと「 $e^{i 2 π / N}$ を $ω_{N}$ に置き換えただけじゃないの？」と思ってしまいがちだが、 1つの係数 $f_{l}$ の $F_{k}$ との間の変換ではなく、 $f_{0}, \dots, f_{N - 1}$ と $F_{0}, \dots, F_{N - 1}$ の間の 線形変換 になっていることが大切。（離散フーリエではなく連続のフーリエ変換も数学的には線形変換）

ベクトルと行列の形で書くと以下になる。

\begin{array}{r} (\begin{array}{c} f_{0} \\ f_{1} \\ f_{2} \\ ⋮ \\ f_{N - 1} \end{array}) = (\begin{array}{ccccc} 1 & 1 & 1 & \dots & 1 \\ 1 & ω_{N} & ω_{N}^{2} & \dots & ω_{N}^{N - 1} \\ 1 & ω_{N}^{2} & ω_{N}^{4} & \dots & ω_{N}^{2 (N - 1)} \\ ⋮ & ⋮ & ⋮ & ⋮ \\ 1 & ω_{N}^{N - 1} & ω_{N}^{2 (N - 1)} & \dots & ω_{N}^{(N - 1) (N - 1)} \end{array}) (\begin{array}{c} F_{0} \\ F_{1} \\ F_{2} \\ ⋮ \\ F_{N - 1} \end{array}) \end{array}

\begin{array}{r} (\begin{array}{c} F_{0} \\ F_{1} \\ F_{2} \\ ⋮ \\ F_{N - 1} \end{array}) = \frac{1}{N} (\begin{array}{ccccc} 1 & 1 & 1 & \dots & 1 \\ 1 & ω_{N}^{- 1} & ω_{N}^{- 2} & \dots & ω_{N}^{- (N - 1)} \\ 1 & ω_{N}^{- 2} & ω_{N}^{- 4} & \dots & ω_{N}^{- 2 (N - 1)} \\ ⋮ & ⋮ & ⋮ & ⋮ \\ 1 & ω_{N}^{- (N - 1)} & ω_{N}^{- 2 (N - 1)} & \dots & ω_{N}^{- (N - 1) (N - 1)} \end{array}) (\begin{array}{c} f_{0} \\ f_{1} \\ f_{2} \\ ⋮ \\ f_{N - 1} \end{array}) \end{array}

離散フーリエ変換を計算するには、長さ $N$ の任意の数列 $a_{0}, a_{1}, a_{2}, \dots, a_{N - 1}$ に対して

b_{k} = \sum_{l = 0}^{N - 1} a_{l} ω_{N}^{k l}, k = 0, 1, 2, \dots, N - 1

となる数列 $b_{0}, b_{1}, b_{2}, \dots, b_{N - 1}$ が計算できればよい

高速フーリエ変換#

$N$ が2の冪乗であれば、離散フーリエ変換の計算が効率化される（具体的には $O (N)$ から $O (N \log N)$ になる）これを使うのが高速フーリエ変換。

単位円での表現#

1の原始 $N$ 乗根 $ω_{N}$ を次々と冪乗すると、複素平面の単位円周上を回転する

Ignoring fixed x limits to fulfill fixed data aspect with adjustable data limits.

../../_images/f5955815a4d112fa0c6ab92173e2812955793fbc3b74c6329f18f8d516d09416.png

定理

$N$ が偶数のとき、次の関係が成り立つ

ω_{N}^{N / 2} = - 1, ω_{N}^{N / 2 + 1} = - ω_{N}, ω_{N}^{N / 2 + 2} = - ω_{N}^{2}, \dots, ω_{N}^{N - 1} = - ω_{N}^{N / 2 - 1}

例えば $N = 8$ のときは、 $ω_{N}^{N / 2} = ω_{N}^{4} = - 1 = - ω_{N}^{0}$ であり、 $ω_{N}^{N / 2 + 1} = ω_{N}^{5} = - ω_{N}^{1} = - ω_{N}$ である。以下同様。

単位円で描くとちょうど反対側に位置する。

そのため、 $N / 2$ 個の $ω_{N}$ についてのみ考えればよい ことがわかる。

Ignoring fixed x limits to fulfill fixed data aspect with adjustable data limits.

../../_images/f494e2caf2bbad8459683cdd04235e0fa40e6edd5f22e786607e46af7454fd3a.png

定理

$N$ が偶数のとき、 $ω_{N}^{2}$ は1の原始 $N / 2$ 乗根として表せる。

ω_{N}^{2} = ω_{N / 2}

前出の定理より、 $ω_{N}^{2 k} = ω_{N / 2}^{k}$ となっている。例えば $ω_{N}^{4} = ω_{N / 2}^{2}$ となっている

Ignoring fixed x limits to fulfill fixed data aspect with adjustable data limits.

../../_images/055bbd7409dd693ee480756cb274167b98f998971fc862f28c66327ee8b1ee73.png

$N - 1$ 次多項式

f (x) = a_{0} + a_{1} x + a_{2} x^{2} + \dots + a_{N - 1} x^{N - 1} = \sum_{l = 0}^{N - 1} a_{l} x^{l}

を定義すると

b_{k} = \sum_{l = 0}^{N - 1} a_{l} ω_{N}^{k l}, k = 0, 1, 2, \dots, N - 1

は

b_{k} = f (ω_{N}^{k}), k = 0, 1, 2, \dots, N - 1

と書くことができる。

つまり、 複素平面上の単位円周の $N$ 等分点で $N - 1$ 次多項式 $f (x)$ を計算すれば $b_{0}, b_{1}, \dots, b_{N}$ が得られる。

この計算を ${FFT}_{N} [f (x)]$ と表す。

{FFT}_{N} [f (x)] = {f (1), f (ω_{N}), f (ω_{N}^{2}), \dots, f (ω_{N}^{N - 1})}

間引き#

上記の多項式は次のように書き直せる

\begin{array}{r} \begin{aligned} f (x) = & a_{0} + a_{1} x + a_{2} x^{2} + \dots + a_{N - 1} x^{N - 1} \\ = & a_{0} + a_{2} x^{2} + a_{4} x^{4} + \dots + a_{N - 2} x^{N - 2} \\ + x (a_{1} + a_{3} x^{2} + a_{5} x^{4} + \dots + a_{N - 1} x^{N - 2}) \\ = & p (x^{2}) + x q (x^{2}) \end{aligned} \end{array}

ただし、

\begin{array}{r} {\begin{cases} p (x) = a_{0} + a_{2} x + a_{4} x^{2} + \dots + a_{N - 2} x^{N / 2 - 1} \\ q (x) = a_{1} + a_{3} x + a_{5} x^{2} + \dots + a_{N - 1} x^{N / 2 - 1} \end{cases} \end{array}

とおいた。この係数の 間引き は高速フーリエ変換における重要な計算手順の一つ。

$p (x^{2})$ は $x$ が複素平面上の単位円周の $N$ 等分点を順にたどって1周するとき $x^{2}$ はそれらを1つおきに進んで2周する 。そのため 前半の $N / 2$ 個のみを計算すればよい から、

{FFT}_{N} [p (x^{2})] = {p (1), p (ω_{N}^{2}), p (ω_{N}^{4}), \dots, p (ω_{N}^{N - 2})}

となる。

また「 $N$ が偶数のとき $ω_{N}^{2} = ω_{N / 2}$ 」という定理より以下のように書き直せる。

{FFT}_{N} [p (x^{2})] = {p (1), p (ω_{N / 2}), p (ω_{N / 2}^{2}), \dots, p (ω_{N / 2}^{N / 2 - 1})}

→ データ数 $N / 2$ 個の離散フーリエ変換 を計算すればいい。

$q (x^{2})$ についても同様に得られる。

{FFT}_{N} [q (x^{2})] = {q (1), q (ω_{N / 2}), q (ω_{N / 2}^{2}), \dots, q (ω_{N / 2}^{N / 2 - 1})}

これら $N / 2$ 個の離散フーリエ変換で得られるものたちを以下のように ${FFT}_{N / 2}$ と表記する

{FFT}_{N} [p (x^{2})] = {FFT}_{N / 2} [p (x)], {FFT}_{N} [q (x^{2})] = {FFT}_{N / 2} [q (x)]

バタフライ#

$f (x)$ は「 $N$ が偶数のとき、 $ω_{N}^{N / 2} = - 1, ω_{N}^{N / 2 + 1} = - ω_{N}, ω_{N}^{N / 2 + 2} = - ω_{N}^{2}, \dots, ω_{N}^{N - 1} = - ω_{N}^{N / 2 - 1}$ 」という定理と $f (x) = p (x^{2}) + x q (x^{2})$ により、以下のように計算できる。

\begin{array}{r} {\begin{cases} f (ω_{N}^{k}) = p (ω_{N / 2}^{k}) + ω_{N}^{k} q (ω_{N / 2}^{k}), & k = 0, 1, \dots, N / 2 - 1 \\ f (ω_{N}^{N / 2 + k}) = p (ω_{N / 2}^{k}) - ω_{N}^{k} q (ω_{N / 2}^{k}), & k = 0, 1, \dots, N / 2 - 1 \end{cases} \end{array}

この計算を バタフライ と呼び、式中の $ω_{N}^{k}$ を 回転因子（ひねり因子） と呼ぶ。

高速フーリエ変換#

${FFT}_{N}$ は係数の間引き（Decimation-In-Time, DIT） $D_{N}$ を行って ${FFT}_{N / 2}$ を実行し、バタフライ $B_{N}$ を施す。

間引きによって ${FFT}_{N / 4}, {FFT}_{N / 8}, {FFT}_{N / 16}, \dots$ と分解していくと最終的に ${FFT}_{1}$ の計算になり、 $0$ 次式（定数）の計算になる。

その後バタフライによって統合していく。

そのため ${FFT}_{N}$ は間引きの繰り返しとバタフライの繰り返しで分割統治法によって計算できる。

この方法を 高速フーリエ変換 （FFT）と呼ぶ。

Pythonの実装#

by chatGPT

import cmath

def fft_recursive(x: list[complex]) -> list[complex]:
    """
    Args: 入力信号（複素数リスト）
    Returns: FFTの結果（周波数成分Fのリスト）（複素数のリスト）
    """
    N = len(x)
    if N <= 1:  # 再帰の終了条件
        return x

    # 偶数番目と奇数番目に分割
    even = fft_recursive(x[0::2])
    odd = fft_recursive(x[1::2])

    # バタフライ操作
    T = [cmath.exp(-2j * cmath.pi * k / N) * odd[k] for k in range(N // 2)]
    return [even[k] + T[k] for k in range(N // 2)] + \
           [even[k] - T[k] for k in range(N // 2)]

FFTの結果:
X[0] = 28.000+0.000j
X[1] = -4.000+9.657j
X[2] = -4.000+4.000j
X[3] = -4.000+1.657j
X[4] = -4.000+0.000j
X[5] = -4.000-1.657j
X[6] = -4.000-4.000j
X[7] = -4.000-9.657j
NumPy FFT結果:
X[0] = 28.000+0.000j
X[1] = -4.000+9.657j
X[2] = -4.000+4.000j
X[3] = -4.000+1.657j
X[4] = -4.000+0.000j
X[5] = -4.000-1.657j
X[6] = -4.000-4.000j
X[7] = -4.000-9.657j

../../_images/94e654589cff875537765955fe23a005febd835fc3bd277945a6bfac845f360f.png

高速フーリエ変換

Contents

高速フーリエ変換#

1の原始N乗根による表現#

高速フーリエ変換#

単位円での表現#

間引き#

バタフライ#

高速フーリエ変換#

Pythonの実装#

1の原始 $N$ 乗根による表現#