PythonでFFT実装!SciPyのフーリエ変換まとめ

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

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する方法(1フレームFFT/DFT)

平均化FFTをする前に、まずは簡単に入力した信号をそのままフーリエ変換する方法を紹介します。

以下にscipyのfftpackを使ったフーリエ変換コードの例を示します。

このコードをそのままコピぺして実行すると、上記の時間波形を周波数分析したグラフが出力されます。

def関数ではspectrumの部分がフーリエスペクトルの結果です。1行で完了してしまうほど簡単ですが、このspectrumは複素数です。

ampで振幅成分、phaseで位相を計算しています。また、freqは周波数軸を作成しています。

以下の図が結果です。上に時間波形、下に周波数波形を示しています。たまにFFTした結果で振幅成分がおかしいネット記事がありますが、振幅成分も生成したsin波の振幅1を示しています。

簡単にFFTした結果

本方法では実はサンプル数はFFTの制約である2の累乗個でなくても機能します(その場合はDFTになるのかな?)

しかし信号の長さが変わると周波数分解能が変わってしまったり、ノイズに弱かったりとなかなか定量評価が難しい所があります。

次は、ある区間毎にフレームを抜き出し、平均化しながらFFTをするためのデータの前処理、そして平均化FFTのコードを紹介します。

平均化FFTのための前処理

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

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

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

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以外のレベルが増加している傾向が明確に表れて来ました。一般的には、この部分をノイズフロアと呼びます。

色々な関数のFFT結果

ここではせっかく覚えたフーリエ変換を使って色々な関数の周波数分析をしてみましょう。

正弦波(\(\sin \theta\))

周波数分析といったらまずは正弦波、というくらいテストされている関数。上記コードのdata部分を以下のように変更する事で今回は振幅1、周波数100Hzの正弦波を作ってみます(ノイズは含ませないパターンです)。

 

100Hzに単一のピークが出てくる結果を得ました。

サイン100Hzの周波数分析結果

ガウシアンノイズ

ガウシアンノイズはガウス分布を持つノイズですが、ホワイトノイズとも呼ばれます。詳しくは「Pythonでガウス分布を持つノイズの作り方と調整方法」で書いていますが、dataを以下のコードに変更するだけです。

以下がガウシアンノイズの周波数分析結果です。正弦波と違って全周波数で特にピークを持たないのが特徴です。

ガウシアンノイズの周波数分析結果

のこぎり波

のこぎり波は「Pythonでのこぎり波を生成!次数の高調波成分を見てみた」で書いていますが、以下のコードで生成可能です。振幅1、周波数100Hzとしています。

以下がのこぎり波の周波数分析結果です。100Hzの基本次数の2倍、3倍…で次数が発生しているのが特徴です。これらの基本次数以外の次数の事を高調波と呼びます。

のこぎり波の周波数分析結果

補足:ミラーリングと周波数軸について

ミラーリングとナイキスト周波数

先ほどまでは周波数分析の結果として、興味帯域しか表示させていませんでしたが、全周波数帯域を表示させると以下の図のようになります(ここではのこぎり波の結果を例に説明します)。

FFT結果のミラーリング(鏡像)ピーク

このように、全周波数帯域を表示させると、みかけ上はサンプリング周波数までデータが存在していますが、サンプリング周波数の半分であるナイキスト周波数\(f_{n}\)を境に折り返されているだけの波形となっている事に気付きます。

このピークの事をミラーリングピークと言ったり、鏡像ピーク、共役鏡像と言ったりするそうです。

そのため、「どの周波数でどれだけの振幅で振動しているか」、といった分析を行う時はこのナイキスト周波数\(f_{n}\)までを扱い、それ以上は切り捨てるという事が行われます。

負の周波数とは?

この記事では便宜上周波数軸をlinspace(0, samplerate, Fs)と正の数値のみで作成していました。

しかし、フーリエ変換は実際は複素数(二乗して負になる数)の計算を行っており、下図のようにフーリエ変換の性質から本来周波数軸は負の値も持ちます(これが本来の姿と思います)。

負の周波数を表現した周波数波形

以下の外部リンクが参考になります(理系なのに数学苦手なので丸投げ!)。

入力データが実数である場合にのみ、FFTの結果が(共役対称のように)ミラーリングされることに注意してください。
厳密に実数の入力データの場合、FFT結果の2つの共役鏡像は、複雑な正弦波の虚数部を打ち消すため、厳密に実数の正弦波(小さな数値の丸めノイズを除く)になります。実際の正弦波。
FFTの結果が共役鏡像化されていない場合、厳密な実数値ではなく、複素数値(非ゼロの虚数成分)を持つ波形を表します。
→詳細は式で。

qstack:FFTが「ミラーリング」されるのはなぜですか?

負の周波数を求めるコード:fftpack.fftfreq()

負の周波数はlinspaceでも作成する事が出来ますが、scipyのfftpackでフーリエ変換しているのであれば、fftpack.fftfreq()を使った方が効率的です。

以下のコードは本記事の一番上に紹介していた単一フレームのフーリエ変換を行う関数ですが、freqの部分をfftpack.fftfreq()を使って書き直しています。

fftfreq()はフーリエスペクトルの長さと信号の時間刻みdを引数として計算されます。

負の周波数が必要になる場面とは?

負の周波数なんて周波数分析で見ないんだからどうでも良いじゃないか!

…と考えがちですが、実は結構重要なんです。

確かに、単純に一番振動している周波数と振幅を分析したいだけであればよくネット記事で語られている正方向の周波数軸がナイキスト周波数分だけあれば問題無いです。

しかし、フーリエスペクトルを使って別途工学演算を行う場合に負の周波数が必要になる事が多々あります。

FFT波形を微積分!Pythonで周波数軸微積分をしてみた」では、フーリエ変換の結果のみを使って微分積分の演算をしています。

この周波数軸演算では角振動数\(\omega\)を使うので、負の周波数と正の周波数を含む正しい周波数軸を計算で使わないと全く異なる演算結果を出してしまいます(この記事を書いた時に僕はそれがわからずハマりました…)。

もしフーリエ変換の結果を使って工学的な計算を行う場合は注意しましょう。

まとめ

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

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

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

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

学んだフーリエ変換を使って様々な関数を周波数分析してみると、それぞれ波形の特徴を捉える事が出来ました。

おまけとして周波数軸に関する注意点を述べました。本ページでPythonのFFTイメージを掴んで頂けたら幸いです。

参考:フーリエ変換の応用例を紹介

逆高速フーリエ変換(IFFT)

フーリエ変換した結果を逆フーリエ変換する事で再び時間波形に戻す事が可能です。

周波数領域でフィルタリングしたり何か処理をした後に時間波形に戻すといった使い道があります。

逆フーリエ変換については「PythonでFFTとIFFT!逆フーリエ変換で時間波形を作る」に詳細を書きましたので参考にして下さい。

2Dフーリエ変換

Pythonで2Dフーリエ変換!画像フィルタリングをする方法」という記事で画像データのフーリエ変換をやってみました(2Dのフーリエ変換)。

本記事の周波数軸に関する説明を読んでからこの記事を読むと、一段と理解の助けになると思います。参考までに。

画像データ(2Dデータ)のFFT結果の例

スペクトログラム分析

スペクトログラムとは時間✖️周波数✖️振幅と3次元的にフーリエ変換を並べたものです。

次数の時間変化を分析する事が出来るため、非常に幅広く使われている分析手法です。是非「Pythonで音のSTFT計算を自作!スペクトログラム表示する方法」をご覧下さい。

スペクトログラムの例

ヒルベルト変換

ヒルベルト変換とはFFTとIFFTを使い時間波形の包絡線を書く事が出来るものです。フーリエ変換を覚えるとこんな事も出来きます。

結構簡単なのでやってみたい方は「Pythonでヒルベルト変換!時間波形の包絡線を求める方法」をご覧下さい。

ヒルベルト変換の例

ケプストラム分析

フーリエ変換の応用として音声分野では声帯情報を分離するケプストラム分析というものがあります。

こちらもフーリエ変換を巧みに使っていますので、ご興味のある方は「Pythonでケプストラム分析!音声から声帯情報を取得してみた」をご覧下さい。

ケプストラム分析の例

周波数応答関数(FRF)とコヒーレンス関数

フーリエ変換を2つの信号に対して行い、クロススペクトルやパワースペクトルを駆使すると入力と出力に関する周波数応答関数(FRF)を算出する事が出来ます。

これも工学的な分析でよく使われるので、コードを確認したい方は是非「Pythonで伝達関数!2つの信号の周波数応答関数を計算する方法」の記事をご覧下さい。

FRFの例

周波数応答関数は同じくフーリエ変換を応用して計算するコヒーレンス関数も併用すると、より確度の高い分析が可能です。

コヒーレンス関数については「Pythonでコヒーレンス関数!FRFのノイズ成分を評価する方法」をご覧下さい。

コヒーレンス関数の例

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

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

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

SNSでもご購読できます。

コメント

  1. アバター スピネル より:

    はじめまして。スピネルと申します。
    現在、フーリエ解析について学んでいる大学生です。
    pythonはまだ初心者でコードを参考にさせて頂いています。
    そこで質問ですが、「import fft_function」の部分でモジュールエラーが起きてしまったので、
    pipコマンドでインストールしようとしたところ、コマンドが見つからないと出てしまいました。
    大変恐縮ですが、「fft_function」について詳しく教えて頂けたら幸いです

    1. wat wat より:

      ご訪問ありがとうございます。
      わかりにくくて申し訳ございません!
      ここでは「fft_function.py」という関数だけを集めたPythonファイルを別途自分で作成し、それをメインのPythonファイルでimportするという方法をとっています。
      そのためpip installは必要ありません。

コメントを残す

*