PythonでFFT!SciPyのFFTまとめ

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

WATLABブログでは時間波形の抽出、オーバーラップ処理、窓関数処理を紹介して来ました。この記事ではそれらのまとめとして、信号処理分野の代表的な計算であるFFT(高速フーリエ変換)を説明します。

こんにちは。wat(@watlablog)です。
いよいよここまで来ました。これまで時間領域の前処理技術を習得して来ましたが、今回はFFTをかけて周波数領域を分析できるようになります

周波数分析概要とFFTのための前処理

時間波形と周波数波形の関係とFFTの威力

当WATLABブログでは信号処理について僕自身が勉強しながら記事を書いています。

信号処理を行う主な目的の一つに周波数分析があります。

周波数とは1秒間に何回振動しているかを数えたものです。そのため、当然時間領域の振動と関係を持っており、図にすると以下のようになります。

単一サイン波の説明図

上図は単一のサイン波\(y(t)=Asin(2\pi ft)\)の時間波形と周波数波形を説明しています。波の周期\(T\)と周波数(振動数)\(f\)には\(f=1/T\)の関係で変換可能であり、振幅\(A\)がわかれば何も難しい理論を要するFFT(Fast Fourier Transform)を使わずとも周波数波形をプロットすることが可能です。

ちなみに周波数波形から時間波形に変換することを逆FFT、またはIFFT(Inverse Fast Fourier Transform)と呼びますが、上図の単一サイン波のように単純なものは紙とえんぴつさえあればプロットできますね。

ではもっと波形を複雑にしてみましょう。複雑な波形でも時間波形と周波数波形の関係は下図のように成立します。

この図の時間波形は、実際にExcelで4つのサイン波を、振幅と位相を変えて足し合わせたものです。

たった4つの波形を合成させただけでも、もはや人間の目と頭では正確に周波数波形を描くことはできません

FFTからは主要なサイン波の個数や、それぞれのサイン波の振幅と位相の情報を得ることができます。これらの情報があると、実際の時間波形がどのような振動成分で構成されているかがわかるので、あらゆる分析に使われています。

例えば、ある構造物の振動原因調査や、どんな材料でできているのかといった工学的な分野に多く使われています。医療現場におけるMRIやCTスキャンの画像処理にもフーリエ変換は活用されています。

今回は以下の図の時間波形を分析してみます。0.5秒程度のデータですが、データサンプル数は12800あります。

波形生成コードは以下です。

波形分割・抽出・オーバーラップ

周波数分析をする場合、計測する対象のデータはある程度一定の振動をしている区間を想定することが多くあります。

実際の振動現象は式で記述するよりも多くのばらつきを含むので、ある程度の時間長分データを計測し、その後に平均化する必要があります。

FFTはそのアルゴリズムからデータ長が2の累乗(\(2^n\))である必要がありますので、平均化に使うデータフレームのサイズも1024や4096…といったものをよく選びます。

この際にオーバーラップという手法を用いる方法と自作したコードを前回記事で紹介しましたので、是非以下の記事を参照下さい。

PythonでFFTをする前にオーバーラップ処理をしよう!

この記事では0.5秒の波形に90%のオーバーラップ率でデータを4096個ずつ抽出した22個の波形データを分析します。以下がオーバーラップ波形です。

オーバーラップ波形

オーバーラップ処理のコード

オーバーラップ処理のコードは以下です。可読性の良さを重視して関数化しています。

窓関数をかける

波形から分析用のデータを抽出したら、窓関数を使ってFFTのための前処理をします。

窓関数については以下の記事で内容とコードを紹介していますので、こちらも必要に応じて参照下さい。

PythonでFFT!SciPyで窓関数をかける

ちなみに、窓関数は自作することも可能です。

Pythonで窓関数が無い場合は?指数窓を自作してみる

先ほどのオーバーラップ波形にハニング窓をかけた結果が以下の図です。

窓関数処理のコード

窓関数処理のコードは以下です。こちらも関数化しています。

平均化FFTのコード

ようやくFFTの出番です。ここでは巷の商用ツールで計算されているようなオーバーラップ処理と窓関数を考慮した周波数平均されたFFTを目指しています。

しかし、このコーディングの最中に筆者自身が直面した壁がありましたので、以下に注意点として残して大きます。まだ真の理解を得たわけではありませんが、今後のキーワード検索のヒントになればと思います。
※理解がクリアになった時点でまた記事にまとめたいと思います。

FFT時の注意点

窓関数を使った場合は補正が必要

窓関数は元の波形を歪ませてしまいます。特に振幅を小さくしてしまうので、FFTで振幅成分を出す時は補正が必要です。

補正については以下の記事に別途まとめましたので、参照下さい。

窓関数使用時の補正!FFTの時に忘れがちな計算とは?

FFTで振幅を出す時に一種の正規化が必要

本ブログは筆者がネットを彷徨い、書籍をあさりながら勉強した結果をアウトプットしています。

FFT計算では主に振幅成分を評価することが多いと思いますが、正しい振幅を得るためにFFTの計算式を「サンプル数の半分で割る」という処理を色々な方が行っていました。コードにすると以下です。

こちらも理論波形を使って検証してみるとやはりこの(\(N/2\))の割り算が無いと正しい振幅成分が得られないという結果を得ました。

その理屈を調べていたのですが、中々僕程度の頭にはピンときませんでした。但し、以下のページの質疑応答がおそらく核心部分と思われますので、引用させて頂きます。

質問

yf = fftpack.fft(y)/(n/2)
で、なぜn/2で割るのかがよくわかりません。

回答
n/2のうちのnはデータ点数で、同じ区間でも解像度が高くなるほど積分値が大きくなってしまうことへの補正です。2は折り返し信号の分をサボったためです。fftshiftを使ったとすれば、正の周波数の負の周波数ではx=0のy軸に対して対称になっていると思います。もう少し細かく見ると、cos項(複素フーリエ係数の実数部分)は正負で同じ値なのに対して、sin項(複素フーリエ係数の虚数部分)は符号が反転していると思います。
少し脱線しましたが、本来ならば同じ周波数の振幅成分が正負に分かれて存在したはずのところを、参考にしている記事では正の部分しか扱っていません。コレ自体はよく行われることです。ただし、本来ならばデータ点数nで割っておけばいいだけのところを、負の周波数成分が無視しているために振幅が合わなくなってしまいます。この点を考慮してnではなくn/2で割っているのです。

teratail:Pythonでフーリエ変換する方法

処理の過程で、サンプル数が多くなると、積分値(離散データなのでΣ計算)が大きくなるというのはイメージが付きます。そのため「Nで割る」という一種の正規化は必要と思いますが、N/2にする意味はフーリエ変換の中身をもう少し理解する必要があります。

ネットに多く溢れているPythonのFFTコードは最後に辻褄を合わせている…ということだと考えられます。高速フーリエ変換の中身に言及する記事が書ければ、再度理解を深めるためにコードを見直したいと思います(できるかな?)。

FFT処理の全コード

前置きが長くなりましたが、以下に今回のFFTを実現する全コードを示します。

ファイルはメイン(main_fft.py)と関数(fft_function.py)に分けて、メインを実行する方法にしました。

まずは関数を集めたファイルです。詳細は各行にコメントを記載しましたので参考にして下さい。

続いてメインの動作を記述したファイルです。後半のグラフ描画域が少し長めですが、基本的に関数を使うことによりメイン計算はシンプルになりました。

このコードを実行すると以下の結果が得られます。

ちなみに、ノイズ成分を0(ランダム波形の係数=0)にすると以下の波形になります。

ノイズ成分を思い切り上げてしまう(ランダム波形の係数=5)と以下の波形になります。このレベルになると、主信号であるサイン波の振幅にノイズが影響を与えてしまうようです。そして基本周波数100Hz以外のレベルが増加している傾向が明確に表れて来ました。一般的には、この部分をノイズフロアと呼びます。

まとめ

このページでは、1フレームしか無い時間波形のFFTをするのではなく、任意長さの時間波形から平均化FFTを求める方法を学びました。

平均化するためには複数のデータフレームが必要で、元の時間波形からフレームをオーバーラップさせながら抽出することで必要な情報をできるだけ含んだ時間波形フレームを作成しました。

窓関数はハニング窓を各時間波形にそれぞれかけ、その時間波形を一つ一つFFT処理しました。このFFT処理時には窓関数による振幅を補正する係数ACFを使いました。

全てのFFTは周波数軸で平均化することにより一つの結果を算出しています。

オーバーラップから窓関数、平均化FFTを自分で実装することができたぞ!
今回の学習で比較的自由に信号を取り扱えるようになってきた気がするので、さらに高度なことにも挑戦しよう!

Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています!

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

SNSでもご購読できます。

コメント

コメントを残す

*