音響専門書籍・文献紹介では、弊社R&D担当者それぞれが日頃参考にしている書籍や、注目している論文などをご紹介します。
紹介した書籍や文献へのリンクも掲載していますので、興味を持たれたらご利用ください。
紹介した書籍や文献へのリンクも掲載していますので、興味を持たれたらご利用ください。
James W. Cooley and John W. Tukey 著
「An Algorithm for the Machine Calculation of Complex Fourier Series」
Mathematics of Computation
19 (90), pp.297-301, 1965
DOI: https://doi.org/10.1090/S0025-5718-1965-0178586-1
今回ご紹介する論文はフリーアクセスとなっています。上記のリンク先ページで、「Full-text PDF」をクリックしていただくと、論文のPDFをダウンロードしていただけます。
キーワード：高速フーリエ変換（FFT）、離散フーリエ変換（DFT）、有限インパルス応答（FIR）
「An Algorithm for the Machine Calculation of Complex Fourier Series」
Mathematics of Computation
19 (90), pp.297-301, 1965
DOI: https://doi.org/10.1090/S0025-5718-1965-0178586-1
今回ご紹介する論文はフリーアクセスとなっています。上記のリンク先ページで、「Full-text PDF」をクリックしていただくと、論文のPDFをダウンロードしていただけます。
キーワード：高速フーリエ変換（FFT）、離散フーリエ変換（DFT）、有限インパルス応答（FIR）
みなさまはFIRフィルタをご存じでしょうか。
FIRとは、Finite Impulse Responseの略で、日本語では有限長インパルス応答と言います。数学的にいえば、あるシステムにインパルス信号を入力したときの出力がインパルス応答で、このインパルス応答を有限の時間で切り取ったものがFIRです。
音響で例えると、お風呂やホールで手をたたいた時（手をたたいた時の音は厳密にはインパルス信号ではありませんが）や、ダミーヘッドマイクの前で手をたたいた時に観測された時間波形がインパルス応答となります。お風呂やホールなどの音響空間や、頭部や外耳などの人体形状というシステムに対して、手をたたいた音というインパルス信号を与えた（入力した）ときの時間波形として観測される出力がインパルス応答ということです。
では、FIRがわかると何がおいしいのでしょうか。
実は、ある線形時不変なシステムにおいては、FIRさえわかれば、任意の入力に対するそのシステムの出力がわかるのです。先ほどの音響で例えるなら、音響空間やダミーヘッドのFIRを利用することで、その音響空間で好きな音を再生（システムに入力）したときに、どのような音波が鼓膜に到来するか（システムからの出力）を計算できます。
なぜそのようなことが可能かというと、システムに入力される音の信号は、たくさんのインパルスが集まったものとみなせるためです。たくさんのインパルスのそれぞれに対して、インパルス応答を計算し、それを重ね合わせることにより、入力された音信号に対する出力を計算することができます。
これは数学的には畳み込みと呼ばれ、次の式で定義されます。 \begin{align*} y[n]=\sum_{k=0}^{N-1} h[k]x[n-k] \end{align*} ここで、$h[t]$は長さ$N$のFIR、$x[n]$と$y[n]$はそれぞれ入力と出力の信号です。$h[k]x[n-k]$はそれぞれのインパルスに対するインパルス応答となり、これを$\sum$によって重ね合わせています。
このように、ある音信号に対してFIRを組み合わせることをFIRフィルタと呼びます。
以上の説明はかなり端折ったものになっていますので、より詳しい部分については専門書等をご参照ください。
FIRフィルタを使えば様々な音響環境での音を再現できるというのは大変魅力的な話ですが、1つ問題点があります。
それは、畳み込みの計算には膨大な計算量を要するという点です。例えば、$44,100\text{Hz}$でサンプリングされた音源とFIRがあり、それぞれの長さは$1$秒だとします。（実際には、FIRは1秒よりも短いことが多いです。）この畳み込みを計算するとき、$44,100$個の量子化された音信号があり、それぞれの音信号$1$個あたりのインパルス応答が$44,100$個ですから、すべてを重ね合わせるには$44,100 \times 44,100$で$1,944,810,000$回の計算が必要になります。 サンプル数が大きくなれば、計算量はこれよりも増えていきます。一般化すると、サンプル数が$N$の音源とFIRの畳み込みに必要な計算量は$\mathcal{O}(N^2)$となります。
いくら現代の計算機（パソコン）が高性能になったとはいえ、このような量の計算を行うのは現実的ではありません。 また近年では、シューティングゲームなど、映像に加えて三次元音響空間で音が到来する方向などを頼りにプレイするゲームが発売されています。これらのゲームでは、三次元音響空間を3Dバイノーラル信号にするための計算が行われており、ヘッドフォンで楽しむことができます。しかしこういったゲームでは、アクションに対する音の遅延を極力少なくすることが求められるため、畳み込みを定義どおりに計算するほどの余裕はありません。
そこで登場するのが、DFT（Discrete Fourier Transform: 離散フーリエ変換）と呼ばれるものです。なぜここでDFTが登場するかというと、それはDFTを用いることで畳み込み計算を簡単に行うことができるからです。DFTは次の式で定義されます。 \begin{align*} A[k]=\frac{1}{N} \sum_{n=0}^{N-1} x[n] W^{-nk}, k=0, 1, \cdots, N-1 \end{align*} ただし、ここでの$x[n]$と$A[k]$は複素数であり、$W=e^\frac{2\pi i}{N}$です。DFTは、時間領域の$N$個のデータ$x[n]$を周波数領域の$N$個のデータ$A[k]$に変換しているとみることができます。 例えるならば、1つの和音の音信号を、どの鍵盤を押したかというデータに変換しているようなイメージです。周波数領域のデータ$A[k]$を時間領域のデータ$x[n]$に戻すには、次の離散フーリエ逆変換を使います。 \begin{align*} x[n]=\sum_{k=0}^{N-1} A[k] W^{nk} \end{align*} さらに、先ほどの畳み込み計算の定義式で出てきた$y[n]$をDFTの定義に沿って計算すると、 \begin{align*} Y[k]&=\frac{1}{N} \sum_{n=0}^{N-1} y[n] W^{-nk}\\ &=\frac{1}{N} \sum_{n=0}^{N-1} \sum_{m=0}^{N-1} h[m]x[n-m] W^{-nk}\\ &=\frac{1}{N} \sum_{m=0}^{N-1} h[m] \sum_{n=0}^{N-1} x[n-m] W^{-nk}\\ &=\frac{1}{N} \sum_{m=0}^{N-1} h[m] \sum_{l=-m}^{N-m-1} x[l] W^{-(l+m)k} (\text{ここで、}l=n-m) \\ &=N\cdot \frac{1}{N} \sum_{m=0}^{N-1} h[m] W^{-mk} \sum_{l=0}^{N-1} x[l] W^{-lk}\\ &=N\cdot H[k]X[k] \end{align*} となります。ただし、$Y[k], H[k], X[k]$はそれぞれ$y[n], h[n], x[n]$のDFTで、また$x[n]$は周期的な信号であるとして、 \begin{align*} x[-l] W^{-(-l)k}=x[N-l] W^{-(N-l)k} \end{align*} を利用しています。つまり循環畳み込みになっています。
この式より、時間領域では$N^2$の計算量が必要だった畳み込み計算が、周波数領域においては$N$回の掛け算で計算できています。
しかしここに落とし穴があります。DFTの定義 \begin{align*} A[k]=\frac{1}{N} \sum_{n=0}^{N-1} x[n] W^{-nk},k=0,1,\cdots,N-1 \end{align*} をよく見ると、ある$A[k]$に対して$N$回の積和演算があります。
$k$は$N$個あるので、全体の計算回数は$N^2$となります。つまり、DFTによって時間領域と周波数領域の変換を行うときの計算量は、もとの畳み込み計算と同じままになっています。このままではDFTを行う意味が無いように思えてしまいます。
前置きが長くなってしまいましたが、ここからが今回の話の本題です。
DFT自体の計算量がもとの畳み込み計算と変わらないという問題を解決するアイデアが、今回ご紹介するCooleyとTukeyによるアルゴリズムです。
アルゴリズム自体はCooleyとTukeyの論文より以前に発見されていたようですが、一般に広めたという功績を鑑みて、彼らの論文についてご紹介します。
彼らのアルゴリズムは次のようになっています。 Complex Fourier seriesを次のように定めます。（彼らの論文に従ってこのように表記します。離散フーリエ変換・逆変換は、定数倍・複素共役の違いを除けばComplex Fourier seriesと同じものとなります。） \begin{align*} x[n]=\sum_{k=0}^{N-1} A[k] W^{nk},n=0,1,\cdots,N-1 \end{align*} これに対して、$N=r_1 \times r_2$を満たす整数$r_1, r_2$が存在するとき、 \begin{align*} n&=n_1 r_1+n_0,n_0=0,1,\cdots ,r_1-1,n_1=0,1,\cdots ,r_2-1\\ k&=k_1 r_2+k_0,k_0=0,1,\cdots ,r_2-1,k_1=0,1,\cdots ,r_1-1 \end{align*} とすることで、もとのComplex Fourier seriesは、 \begin{align*} x[n]=x[n_1, n_0]&=\sum_{k=0}^{N-1} A[k] W^{nk}\\ &=\sum_{k_0=0}^{r_2-1} \sum_{k_1=0}^{r_1-1} A[k_1 ,k_0] W^{nk_1 r_2} W^{nk_0} \end{align*} と変形できます。さらに、オイラーの公式より、 \begin{align*} W^{nk_1 r_2}=W^{n_0 k_1 r_2} \end{align*} を使うことで、 \begin{align*} A_1 [n_0, k_0] &=\sum_{k_1=0}^{r_1-1} A[k_1 ,k_0] W^{n_0 k_1 r_2}\\ x[n_1, n_0] &= \sum_{k_0=0}^{r_2-1} A_1 [n_0 ,k_0] W^{(n_1 r_1+n_0)k_0} \end{align*} と2段階に分割できます。
1段階目では$[n_0, k_0]$の$N$個の組み合わせに対して$r_1$個の総和をとっているため、（$W$の累乗の計算を別にすれば）計算回数は$N r_1$となり、同様にして2段階目の計算回数は$N r_2$となるため、合計の計算回数は$N(r_1+r_2)$となります。
さらに、$N=r_1 \times r_2 \times \cdots \times r_m$と分解できるとき、Complex Fourier seriesの計算を再帰的に分割することで、合計の計算回数は$N(r_1+r_2+\cdots +r_m)$となり、$N=r^m$のときには計算回数は$Nr \log_r N$となります。
このように、計算を分割することで計算量をもともとの$N^2$より減らすことができるというところに感動したので、今回このアルゴリズムを説明した論文を紹介させていただきました。
このアルゴリズムですごいと感じたところは計算量の減少だけではありません。次の3つの点にも感動しました。
まず、各段階の中では、総和をとるパラメータ以外は独立であるため、並列に計算できること。
次に、計算中にデータを保存するために必要な配列の長さは$N$であること。（この論文が書かれた当時はトランジスタを使った計算機が登場して間もないような時代だったため、メモリの量も今とは比べ物にならないほどの小ささでした。そのため、メモリの無駄がないというのはとても大きなメリットでした。）
そして、$N=2^m$のとき、データを$m$ビットのインデックスで指定することができ、さらに$l$段階目$(1 \leq l \leq m)$の計算には、$(m-l)$ビット目のみが異なる2つのインデックスのデータを使えばよいため、データインデックスの指定が容易なことです。個人的には、この点が一番美しいと思っています。
以上が論文の要旨です。これこそが、DFTを高速で計算する手法である、Cooley-Tukey FFT アルゴリズム（FFT: Fast Fourier Transform、高速フーリエ変換）です。今までに、実数値データに対するFFTや、分散メモリ型計算機向けのFFTなど、多くのFFTアルゴリズムが開発されてきましたが、そのほとんどはCooley-Tukey FFTアルゴリズムがもとになっています。
現在使われているFFTライブラリの多くは、データ数$N$を$N=2^m$としているものが多く、そのときの計算量は$\mathcal{O}(N \log_2 N)$となります。例えば入力データ数が$N=2^{15}=32,768$であれば、$N^2$は$1,073,741,824$となるのに対し、$N log_2 N$は$491,520$となりますから、計算量が劇的に改善されていることがわかります。
冒頭で説明したFIRフィルタの畳み込み計算も、FFTを使うことで$\mathcal{O}(N \log_2 N)$の計算量となり、とても実用的なものとなります。現在では、FIRフィルタはデジタルオーディオ機器やコンピューターゲームなど様々な機器やアプリケーションで利用されています。
我々が普段利用しているアプリケーションの向こう側には、美しい理論が隠れています。
今回はそんな中の1つ、CooleyとTukeyのアルゴリズムについてのお話でした。
η
