Pythonで高速フーリエ変換をCooley-Tukey法で書く

  • このエントリーをはてなブックマークに追加

PythonにはNumPyやSciPyによる優秀な高速フーリエ変換(FFT)アルゴリズムがありますが、基本的なFFTの仕組みを学習するにもPythonは最適です。ここでは最も一般的なFFTアルゴリズムであるCooley-Tukey法をPythonのNumPyで記述し、numpy.fftの結果と比較します。さらに、データ点が2の冪乗でない場合のゼロパディング処理も実装します。

こんにちは。wat(@watlablog)です。ここではCooley-Tukey法による高速フーリエ変換をPythonで書き、そのアルゴリズムを勉強します

DFTとFFTの違い

DFT:離散フーリエ変換

WATLABブログでは以下の記事で離散フーリエ変換(DFTDiscrete Fourier Transform)を紹介しました。
Pythonで検証しながらSwiftで離散フーリエ変換を実装する

この記事では基本的なDFTの式とPythonやSwiftによる実装例を紹介しましたが、「実用的ではない」と結論付けています。その理由はDFTの計算量が\(O(N^{2})\)とデータ数の二乗に比例して膨大な計算が必要になってしまうことにあります。

プログラム的にはとても簡単に書けてしまいますが、音声データ処理といったkHzオーダーのサンプリングレートを持つデータの処理をしようとすると一般のPCではかなり長い時間待たされてしまうことでしょう。

その計算量に関する問題を改善するために、高速フーリエ変換が考案されました。

FFT:高速フーリエ変換

FFTの特徴

高速フーリエ変換(FFTFast Fourier Transform)とは、離散フーリエ変換(DFT)を高速に計算するために考案された手法です。

線形かつ周期的なデータである必要はありますが、その計算量は\(O(N \log N)\)になるためデータが膨大になるに従いDFTと比べものにならないくらいの計算速度改善を実感します。

Cooley-Tukey法

FFTの一般的なアルゴリズムにCooley-Tukey法があります。この手法は問題をより小さな部分問題に分解する「分割統治」戦略をとっています。

アルゴリズムは①データの分割②再帰的なDFT③データの結合バタフライ演算)の順番で計算を行うのが特徴です。

数式による説明は英語版のWikipedia\(^{[1]}\)がわかりやすいと思います。
筆者はそんなに詳しく説明できないですが、要は重複計算をしてしまうDFTを最適化して効率良く必要な計算のみしていこうという方針のようです。単純に計算数を少なくすることで高速化しています。

この記事では代表的なFFTアルゴリズムであるCooley-Tukey法をPythonでコーディングして結果を確認します。ただし筆者も勉強中であり、ChatGPT-4に数学やコーディングを教えてもらいながら学習的な用途で書いたものなので、この記事のコードは計算高速性について最適化されたものではありません

FFTは高速性を求めると非常に奥が深い計算です。必要であれば他の研究者のコードや安定したライブラリを使いましょう。

ちなみにPythonライブラリを多用したコードは以下の記事をご覧ください。
PythonでFFT実装!SciPyのフーリエ変換まとめ
ただPythonでcsvから離散フーリエ変換をするだけのコード

動作環境

この記事のコードは以下の環境で実行しました。

PC

■Mac
                                          
Mac OSmacOS Ventura 13.2.1
CPU1.4[GHz]
メモリ8[GB]

Python

Python Python 3.9.6
PyCharm (IDE) PyCharm CE 2020.1
Numpy 1.21.1
matplotlib 3.4.3

Cooley-Tukey法によるFFTコード

データ長が2の冪乗の場合

関数部のPythonコード

以下にCooley-Tukey法で書いたFFTのPython関数を示します。標本点数が1になるまで再帰的に同じ関数を使い続ける構造となっていることがわかると思います。

再帰的な処理はプログラミングでよく出てきますが、WATLABブログでは基本情報技術者試験の勉強をした時にマージソートの記事で説明しましたので参考にしてみてください。

この関数はデータ長が2の冪乗(2, 4, 6, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192, 16384, 32768…)である場合のみ使用できます。これはデータをevenとoddに再帰的な分割をしていることから来る制約です。

結果評価までのサンプルコード(コピペ動作OK)

下のコードはコピペで動作するサンプルコードです。sample_wave()関数でデータポイントを指定した直流成分付き正弦波を生成し、自作のDFT自作のFFT(上記関数)、numpy.fftを計算して計算時間を確認します。自作のFFTは再帰的処理をしている都合、振幅計算はfft_post()で行い、その結果をmatplotlibで表示させます。

実行結果:計算速度とnumpy.fftとの比較検証

プログラムを実行するとコンソールには各手法の計算時間が表示されます。今回は512データポイントの場合ですが、自作DFTが0.55[s]、自作FFTが0.0098[s]、numpy.fftが0.00095[s]という結果でした。
numpyが圧倒的に高速なのはさすがですが、自作FFTも自作DFTより数十倍も速くなっています

コードの最後で「np.allclose(F, F_numpy)」を使って自作FFTとnumpy.fftの結果を比較しています。np.allcloseは全ての要素が一致している場合Trueを返しますが、先ほどの結果もTrueになっており自作したFFTとnumpy.fftの結果が一致していることを確認しました。

実行結果:時間波形と周波数波形

下図がグラフです。周波数波形には鏡像ピークが出ていますが、これが正常です。通常振幅評価をする時はナイキスト周波数で表示を切るのが一般的ですが、今回は全体を俯瞰して検証する目的でそのまま表示させています。

fft_post()でスケーリング係数sfを分岐してかけていますが、この操作のおかげで直流成分も指定した数値となっていることも確認できました。この辺の話は以下の記事に書いた通りです。
Pythonで検証しながらSwiftで離散フーリエ変換を実装する

データサンプル数と計算時間の関係

DFTとFFTを比較するなら、やはりデータサンプル数と計算時間の関係をプロットしたくなります!こちらがその結果です。

音声データの周波数分析はフレームサイズが2048より大きく、さらにオーバーラップで何回も計算することを考えるとFFTを使わない手はないというほどの差ですね…。

DFTとFFTのデータサンプル数と計算時間の関係図

自作FFTもかなり良い…と思いきや、numpy.fftはさすがの計算速度!素直にnumpyやscipyを使いましょう。

Pythonでは素直にnumpyやscipyを使いましょう(もう一回言っておきます)。
ただし、ライブラリが使えない状況(Excel VBAやその他言語)の場合は一度フルスクラッチしておけばいつでも使える状態になるので、このようなコードの写経はしておいた方が良いと思いました!

データ長が2の冪乗でない場合:ゼロパディング法

ゼロパディングとは?

Cooley-Tukey法ではデータ長が2の冪乗数である必要がありましたが、その制約を取ることも可能です。
ゼロパディングを行うことでそれを実現しますが、ゼロパディングについては過去WATLABブログでも周波数分解能を操作する時に使ったことがあります。詳細は以下の記事をご確認ください。
ゼロパディングFFTで高周波数分解能にするPythonコード

今回行うゼロパディングとはデータ長を2の冪に合わせるために0配列を水増しすることです。ゼロを追加するのでフーリエ変換した時の振幅スケーリング時に辻褄を合わせる必要があります。

ゼロパディング付きCooley-Tukey法のPythonコード(コピペ動作OK)

こちらがゼロパディング付きCooley-Tukey法のFFTコードです。とは言え自作なのであくまで参考までにどうぞ。ゼロパディングによる振幅補正係数acfは先ほど紹介した記事に記載のまま使っています。

実行結果:時間波形と周波数波形

データポイント数Nを511と外していますが、エラーにならない結果を得ました。511点がゼロパディングにより512点になりました。

511点データに対するFFT結果例

1800点のデータでもトライ。2048点にゼロパディングされてエラーなしで終了です。

1800点データに対するFFT結果例

振幅補正係数を使っているのでそれなりに正しい振幅が得られている傾向はありますが、ややゼロパディングすると小さい振幅になっているような…。ここはもしかしたらバグがどこかにありそうです(まだ探せていないのでToDoにしておきます)。

データサンプル数と計算時間の関係

2の冪乗データ数をわずかに超えると、ゼロパディングにより次の冪乗数になるため容易に理解できる結果となりました。

ゼロパディング付きFFTのデータサンプル数と計算時間の関係

調査用コードも貼っておきます。

まとめ

今まで「高速フーリエ変換は難しくてライブラリを使わないで書くのは無理」と食わず嫌いしていましたが、Cooley-Tukey法であればそんなにハードなプログラミングを要求されないことがわかりました。

この記事ではDFTとFFTの違いを解説し、有名なCooley-Tukey法によるFFTコードを紹介、DFTと比べて圧倒的に計算時間が速くなることも確認しました。

また、Cooley-Tukey法ではデータサンプル数が2の冪乗でないといけませんが、ゼロパディングを使うことでその制約を排除できることを確認。これでその他の言語が信号処理ライブラリを持っていない場合でも高速フーリエ変換を実装することができそうです。
(そしてPythonではnumpyやscipyを使った方が良いというのは改めて実感しました)

参考文献

[1]:Wikipedia(E), Cooley–Tukey FFT algorithm

高速フーリエ変換(FFT)を再学習しました!
Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています!

SNSでもご購読できます。

コメントを残す

*