logo

高速フーリエ変換アルゴリズム 📂フーリエ解析

高速フーリエ変換アルゴリズム

概要1

離散フーリエ変換DFTは、数式的な定義に従って計算すると、O(N2)\mathcal{O}(N^{2})の時間計算量を持ちますが、以下で説明するアルゴリズムに従って計算すると、時間計算量がO(Nlog2N)\mathcal{O}(N\log_{2}N)に低減します。この高速フーリエ変換を使用して離散フーリエ変換を高速に実行することができます。これは高速フーリエ変換fast Fourier transform, FFTと呼ばれています。

構築

2つの数字を掛けて、それを別の数字に加える操作を11回の演算operationと呼びましょう。すると、i=0n1anbn\sum\limits_{i=0}^{n-1}a_{n}b_{n}の値を求めるにはnn回の演算が必要です。

n=00anbb=a0b0=0+a0×b01 operationn=01anbb=a0b0+a1b1=(0+a0×b01 operation)+a1×b12 operationsn=02anbb=a0b0+a1b1+a2b2=((0+a0×b01 operation)+a1×b12 operations)+a2×b23 operations \begin{align*} \sum\limits_{n=0}^{0} a_{n}b_{b} &= a_{0}b_{0} = \overbrace{0 {\color{red} +} a_{0} {\color{red} \times} b_{0}}^{\color{red} 1 \text{ operation}} \\ \sum\limits_{n=0}^{1} a_{n}b_{b} &= a_{0}b_{0} + a_{1}b_{1} = \overbrace{\big( \overbrace{0 {\color{red} +} a_{0} {\color{red} \times} b_{0}}^{\color{red} 1 \text{ operation}} \big) {\color{#5882FA} + } a_{1} {\color{#5882FA} \times} b_{1}}^{\color{#5882FA}2 \text{ operations}} \\ \sum\limits_{n=0}^{2} a_{n}b_{b} &= a_{0}b_{0} + a_{1}b_{1} + a_{2}b_{2} = \overbrace{\bigg( \overbrace{\big( \overbrace{0 {\color{red} +} a_{0} {\color{red} \times} b_{0}}^{\color{red} 1 \text{ operation}} \big) {\color{#5882FA} + } a_{1} {\color{#5882FA} \times} b_{1}}^{\color{#5882FA}2 \text{ operations}} \bigg) {\color{#FE9A2E} + } a_{2} {\color{#FE9A2E} \times} b_{2}}^{\color{#FE9A2E}3 \text{ operations}} \\ \end{align*}

さて、離散フーリエ変換の定義を思い出してみましょう。

線型変換 FN:CNCN\mathcal{F}_{N} : \mathbb{C}^{N} \to \mathbb{C}^{N}離散フーリエ変換と呼びます。

FN(a)=a^=[a^0a^1a^N1],a^m=n=0N1ei2πmn/Nan(0m<N)(1) \mathcal{F}_{N}(\mathbf{a}) = \hat{\mathbf{a}} = \begin{bmatrix} \hat{a}_{0} \\ \hat{a}_{1} \\ \dots \\ \hat{a}_{N-1} \end{bmatrix} ,\quad \hat{a}_{m} = \sum_{n=0}^{N-1}e^{-i2\pi mn /N} a_{n}\quad (0\le m < N) \tag{1}

このとき、a=[a0a1aN1]T\mathbf{a} = \begin{bmatrix} a_{0}& a_{1}& \dots& a_{N-1} \end{bmatrix}^{T}です。

a^m\hat{a}_{m}を計算するにはNN回の演算が必要で、a^\hat{\mathbf{a}}を計算するにはこれをNN回繰り返す必要があるため、離散フーリエ変換を計算するには合計でN2N^{2}回の演算が必要です。つまり、O(N2)\mathcal{O}(N^{2})時間計算量を持っています。これは、コンピュータ計算の観点からフーリエ変換がかなりのコストを要することを意味します。

アルゴリズム

データの長さNN合成数N=N1N2N = N_{1}N_{2}としましょう。そして、インデックスm,nm, nを次のように定義します。

m=mN1+m,n=nN2+n m = m^{\prime}N_{1} + m^{\prime \prime},\quad n = n^{\prime}N_{2} + n^{\prime \prime}

すると、0m,nN210 \le m^{\prime}, n^{\prime \prime} \le N_{2}-1および0m,nN110 \le m^{\prime \prime}, n^{\prime} \le N_{1}-1です。(1)(1)の指数部分を次のように表現できます。

ei2πmn/N=ei2π(mN1+m)(nN2+n)/(N1N2)=ei2π(mnN1N2+mnN1+mnN2+mn)/(N1N2)=ei2πmnei2π((mn/N2)+(mn/N1)+(mn/N))=ei2π((mn/N2)+(mn/N1)+(mn/N)) \begin{align*} e^{-i2\pi mn /N} &= e^{-i2\pi (m^{\prime}N_{1} + m^{\prime \prime})(n^{\prime}N_{2} + n^{\prime \prime})/(N_{1}N_{2})} \\ &= e^{-i2\pi (m^{\prime}n^{\prime}N_{1}N_{2} + m^{\prime}n^{\prime \prime}N_{1} + m^{\prime \prime}n^{\prime}N_{2} + m^{\prime \prime}n^{\prime \prime})/(N_{1}N_{2})} \\ &= e^{-i2\pi m^{\prime}n^{\prime}} e^{-i2\pi \big( (m^{\prime}n^{\prime \prime}/N_{2}) + (m^{\prime \prime}n^{\prime}/N_{1}) + (m^{\prime \prime}n^{\prime \prime}/N) \big)} \\ &= e^{-i2\pi \big( (m^{\prime}n^{\prime \prime}/N_{2}) + (m^{\prime \prime}n^{\prime}/N_{1}) + (m^{\prime \prime}n^{\prime \prime}/N) \big)} \\ \end{align*}

これを(1)(1)に代入すると、

a^m=n=0N1ei2πmn/Nan=n=0N21n=0N11ei2π((mn/N2)+(mn/N1)+(mn/N))anN2+n=n=0N21[n=0N11ei2πmn/N1anN2+n]ei2π[(mn/N2)+(mn/N)] \begin{align*} \hat{a}_{m} &= \sum_{n=0}^{N-1}e^{-i2\pi mn /N} a_{n} \\ &= \sum_{n^{\prime \prime}=0}^{N_{2}-1} \sum_{n^{\prime}=0}^{N_{1}-1} e^{-i2\pi \big( (m^{\prime}n^{\prime \prime}/N_{2}) + (m^{\prime \prime}n^{\prime}/N_{1}) + (m^{\prime \prime}n^{\prime \prime}/N) \big)}a_{n^{\prime}N_{2}+n^{\prime \prime}} \\ &= \sum_{n^{\prime \prime}=0}^{N_{2}-1} \left[ \sum_{n^{\prime}=0}^{N_{1}-1} e^{-i2\pi m^{\prime \prime}n^{\prime}/N_{1}}a_{n^{\prime}N_{2}+n^{\prime \prime}} \right] e^{-i2\pi [ (m^{\prime}n^{\prime \prime}/N_{2}) + (m^{\prime \prime}n^{\prime \prime}/N) ] } \end{align*}

上記の式に従うと、各括弧内の[n=0N11]\left[ \sum_{n^{\prime}=0}^{N_{1}-1} \right]を計算するのにN1N_{1}回の演算、括弧の外のn=0N21\sum_{n^{\prime \prime}=0}^{N_{2}-1}を計算するのにN2N_{2}回の演算が必要です。したがって、a^m\hat{a}_{m}を計算するには合計で(N1+N2)(N_{1} + N_{2})回の演算が必要です。a^\hat{\mathbf{a}}を得るにはこれをNN回繰り返す必要があるため、合計でN(N1+N2)N(N_{1} + N_{2})のコストがかかり、N2N^{2}よりも減少することが確認できます。

括弧内を注意深く見ると、N1N_{1}が再び合成数の場合、同じロジックを適用できることがわかるでしょう。したがって、データの長さがN=N1N2NkN = N_{1} N_{2} \cdots N_{k}といった合成数の場合、時間計算量は次のように減少します。

O(N2)O(N(N1+N2++Nk)) \mathcal{O}(N^{2}) \searrow \mathcal{O}\big( N(N_{1} + N_{2} + \cdots + N_{k}) \big)

ここで、NN22のべき乗N=2kN = 2^{k}と仮定してみましょう。すると、log2N=k\log_{2}N = kであり、N2=2kN^{2} = 2^{k}から2k(2k)2^{k}(2k)だけ減少するため、時間計算量は次のように減少します。

O(N2)O(2Nlog2N) \mathcal{O}(N^{2}) \searrow \mathcal{O}(2N \log_{2}N)

補足

これは1965年にCooleyとTukey2によって提案されたため、Cooley-Tukeyアルゴリズムとも呼ばれています。ただし、彼らが最初に発明したわけではありません。ガウスも同様のアルゴリズムを研究しましたが、正しく発表しなかったため、この事実は後に明らかになりました3


  1. Gerald B. Folland, Fourier Analysis and Its Applications (1992), p252-253 ↩︎

  2. J. W. Cooley and J. W. Tukey, An algorithm for the machine calculation of complex Fourier series, Mathematics of Computation 19 (1965), 297-301. ↩︎

  3. M. T. Heideman, D. H. Johnson, and C. S. Burms, Gauss and the history of the fast Fourier transform, Archive for the History of the Exact Sciences 34 (1985), 264-277. ↩︎