はじめに
鈍行ではありますが、高速フーリエ変換の修得を進めています。先日、ChatGPTを利用しつつFFTするRustの実装ができました。これは、2021年9月号Interface 星野秋人さんの記事と「道具としてのフーリエ解析」を参考にしています。改めて、ここまでの理解で高速フーリエ変換について記述します。
高速フーリエ変換とは
高速フーリエ変換(fast Fourier transform, FFT)は、離散フーリエ変換(O(N^2))を効率的に計算するアルゴリズム(O(NlogN))です。この高速フーリエ変換により、コンピュータに現実的な処理速度で実装可能になります。ただし、N(入力データ数)が2の累乗である必要があります。これは高速フーリエ変換が処理対象データを再帰的に1/2個ずつに分割するアルゴリズムであることに由来します。
離散フーリエ変換については、こちらで説明しています。
高速フーリエ変換のアルゴリズム
4Stepで記述したアルゴリズムが以下になります。
- Step1 バタフライ演算を行う
- Step2 データ列(N個)をN/2個に分割する
- Step3 データ数が2個になるまで、Step1とStep2を繰り返す
- Step4 ビットリバースを行う
なお、バタフライ演算とビットリバースは高速フーリエ変換の中核をなす制御です*1。
Step1 バタフライ演算を行う
まず、バタフライ演算とは、2つのデータに2種類の演算をし、2つのデータを求める演算です。例えば、(a,b)という入力に対して、(a+b)と(a-b)という2値出力を得るなどです。この演算過程の可視化としてシグナルフロー図というものがあります。「シグナルフロー図 FFT」で検索するといろいろ出てきます。
FFTのバタフライ演算は、2つのデータ(a,b)から「加算したデータ(a+b)」と「減算し、回転子Wのk乗をかけたデータ(W^k * (a-b))」の2値を求める演算を行います。2つのデータの選び方は、データ列を前半と後半に二分割して、それぞれの先頭から組にしていきます。kは組のインデックスになります。また、回転子Wの詳細は、高速フーリエ変換の前段である離散フーリエ変換を理解するを参照ください。なお、基底(データが2個だけの場合)は、「加算したデータ(a+b)」と「減算したデータ(a-b)」の2値を求める演算となります。回転子Wが1または-1になるためです。
演算結果は、「加算したデータ(a+b)」を前半分に、「減算し、回転子のk乗をかけたデータ(W^k * (a-b))」を後半分に、格納していきます。
Step2 データ列(N個)をN/2個に分割する
Step1の演算により、データ列の前半が加算データ群、後半が回転子*減算データ群に分割されました。この前後半に二分割したデータ列を用意します。
Step3 データ数が2個になるまで、Step1とStep2を繰り返す
データ数が2個になるまで、Step1とStep2を再帰処理します。Step2により、データ数が半分になっていくので、いつかはStep2でのデータ列が2個になります。この二分割再帰処理でフーリエ変換演算を完了させるため、離散フーリエ変換に比べて計算量を減らすことができます。
Step4 ビットリバースを行う
ビットリバースにより出力データを並び替えて0,1,2,...,N-1の順番に戻します。ビットリバースとは、数を2進数表示しそれを逆順にしたときの値、です。例えば、3bitのビットリバースは以下になります。
| 値 | ビットリバース | 2進数 | 逆転 |
|---|---|---|---|
| 0 | 0 | 000 | 000 |
| 1 | 4 | 001 | 100 |
| 2 | 2 | 010 | 010 |
| 3 | 6 | 011 | 110 |
| 4 | 1 | 100 | 001 |
| 5 | 5 | 101 | 101 |
| 6 | 3 | 110 | 011 |
| 7 | 7 | 111 | 111 |
FFTでは、バタフライ演算により入力データの各要素ごとにDFTを行った結果が得られます。しかし、その結果の順番が0,1,2,3,...,(N-1)とはなっておらず、入れ替わりが発生しています。そのため、バタフライ演算の結果順番を0,1,2,3,...,(N-1)に並び替える必要があります。入れ替わりが発生していますが、規則性があり、ビットリバース順に並べ替えることで期待順番0,1,2,3,...,(N-1)に戻せます。
以上で高速フーリエ変換完了です。
高速フーリエ変換の実装
github.com
にてRustで実装しています。
なぜ高速化できるのか？
ここまでで、実装に必要なアルゴリズム説明をしました。しかし、離散フーリエ変換から高速フーリエ変換が実現できる理由を説明していません。私は良い感じに説明をできるほどの理解はできておりません。ざっくり私の理解だと、離散フーリエ変換のDFT行列の要素は回転子Wのk乗で構成され、同じ計算をしている箇所があり、それらを"まとめる"ことで、高速化を実現している、"まとめる"方法はバタフライ演算に対応する、という感じです。詳細は参考文献を参照ください。
おわりに
高速フーリエ変換の概要とアルゴリズムをまとめました。高速フーリエ変換アルゴリズムの中身にバタフライ演算とビットリバースというものがあり、さらに再帰処理をしている、ということを今回初めて知った気がします。ようやくここまで知ることができました。今後は、波形データへのFFT、マイコンでFFTあたりをしたいなと思います。あとAtCoderの高速フーリエ変換も理解したいです。まだまだ高速フーリエ変換と戯れる感じになりそうです。
*1:これらについて「道具としてのフーリエ解析」に詳述されており、長い時間をかけて理解に至れました。
