PythonでFFT波形から任意個数のピークを自動検出する方法

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

科学技術計算では波形からピークを読み取る必要性に迫られます。Pythonではscipyのsignal.argrelmaxを使って簡単にピーク検出が可能です。ここではFFT波形を例に自動でピーク検出する方法を紹介します。

こんにちは。wat(@watlablog)です。ここではFFT波形からピークを自動検出する方法を紹介します

ピーク検出の勘所

FFT波形はピークの周波数・振幅情報が重要

フーリエ変換して得られる周波数波形は、周期的な振動成分の情報を明らかにします。

以下の図に周波数波形の例を示します。このようにある周波数で振動している現象を計測した場合、ピーク(Peak)が観測されます。

このピークには、どれくらいの振動数(周波数成分)、大きさ(振幅成分)で振動しているかといった情報を含んでおり、計測結果から現象を推定するのに役立てることができます。

FFTのピーク

大抵の場合、ピークは1つだけではありません。ピーク検出アルゴリズムがあれば下図のように複数のピークを一度に検出することができ、振幅成分と周波数成分を自動記録していくことが出来ます。

複数のピーク

一度だけの実験や解析であれば、手作業で記録しても良いですが、回数が増えてくると大変な手間です!大量作業はプログラムに任せましょう!

ピーク検出はノイズ対策が必要

ピークを検出するには課題もあります。多くの実験データは外乱によるノイズを含みます。

ノイズを含むデータをフーリエ変換した場合、下図のように本来検出したいピークとは別に誤検出してしまうピークもあります。

ノイズを含み誤検出したピーク

感度良くピークを検出してしまうとこのように誤検出ピークが多くなり、これでは情報の整理がままなりません。

結局手作業で選別するならピークを検出した意味がない!

そのため、通常はあるをもってピークを検出することで感度を鈍らせるノイズ対策処理を行います(下図はイメージ図)。

ピーク検出は元波形の微分波形が0になる(元波形の傾きが0)点を探すというアルゴリズムを基本にしており、上記ノイズ対策のために多項式近似を使ったり、上側や下側のピークのみ検出するために方向を特定する処理が入っていたりと、やや複雑です。

汎用的なピーク検出関数を自作するのは出来なくもないですが(…多分)、Pythonには自分で微分とかしなくても良いピーク検出のメソッドが既にあるため、このページではその方法を紹介するのみとします。

scipy.signal.argrelmaxを使う

scipyのsignalに入っているargrelmaxを使えば簡単にピーク検出が可能です。

新しい外部ライブラリを利用する時は、まず最初に1次情報として公式ドキュメントを確認する事をオススメします。そのため公式ドキュメントへのリンクを以下に示します。

公式ドキュメント:scipy.signal.argrelmax

それでは以下より実際にPythonのscipy.signal.argrelmaxを使ってピーク検出するコードの例を紹介していきます。

Pythonでピーク検出するコード

今回はサンプルの時間波形をフーリエ変換し、周波数波形から任意個数のピーク情報を抽出するプログラムを書いていきます。

import文

importは配列処理、演算でおなじみのnumpy、ピーク検出とフーリエ変換用のscipy、グラフで確認するためのmatplotlibを使います。

フーリエ変換する関数

サンプル波形を以下のコードで作成します。サンプル波形はsignal.sawtoothで基本周波数1[Hz]ののこぎり波(高調波が綺麗に出るため)を作成し、その波形にnp.random.normalで生成したホワイトノイズを重畳させています。
のこぎり波については「Pythonでのこぎり波を生成!次数の高調波成分を見てみた」をご覧下さい。

ここではフーリエ変換については詳しく語りませんが、気になる方は「信号処理」カテゴリのフーリエ変換の項目をご確認下さい。

ピーク検出する関数

以下が本題のピーク検出関数です。findpeaks(MATLABと同じにしてみました)という名前のdef関数で作成しており、入力波形の横軸x、縦軸の値y、抽出するピーク数n、ピーク感度に関する幅(orderとも呼ぶ)wを引数とし、returnとしてピークの位置index、ピークの縦軸値peaksを返します。

この関数は周波数波形だけでなく、xとyに時間軸と値を渡せば時間波形でも機能します。

関数は以下のコードで実行します。

グラフ表示までの全コード(コピペ用)

コピペ用に全コードを載せておきます。ここではグラフ表示まで載せているので、このコードを丸ごとコピペすれば動作すると思います。

実行結果

実行結果を説明していきます。ここではfindpeaks関数に渡すw値による違いを見ていきます。

w=2

まず全コードに載せているようにw=2の結果ですが、以下の図のように、周波数波形(下図)では基本周波数とその倍数成分(高調波)以外にも多くのピークが誤検出されていることがわかりました。

ピーク検出結果(w=2)

w=6

次にw=6の結果です。この結果は基本周波数とその高調波成分のみがピークとして検出されています。これで目的達成です。

ピーク検出結果(w=6)

おまけ:負のピークを検出する方法

実はargrelminを使うと負のピーク(下側に凸のピーク)を検出することが可能です。

以下のdef関数はdirectionという引数を使ってargrelmaxを使った正のピーク検出とargrelminを使った負のピーク検出を選択できるように元の関数を改良したものです。

以下のコードのようにdirectionに文字列として「"+"」を渡せばargrelmaxによる正のピーク検出、「"-"」を渡せばargrelminによる負のピーク検出を行うようになります。

以下が実行結果です。サンプル波形は負のピークがノイズしかありませんが、しっかり下側に凸のピークを検出していることが確認できました。

負のピーク検出結果(w=6)

さらにノイズを軽減するためには?

ノイズをどんな状況でも取り切るw(argrelmax/argrelminではorder)の値はありません。実際のデータで最適な値を決めるべきハイパーパラメータです。

実際の計測ではノイズの度合いも毎回変わってしまうこともよくあることで、以外と自動化する時にはノイズ対策が難しくなります。

しかし、状況によってはさらにノイズを軽減させることも十分可能です。

計測対象が回転振動であれば、回転数のデータを使ってピーク探索範囲を限りなく狭めることができたり(ピーク検出アルゴリズムが必要無いくらいに)、ノイズの原因が高周波成分や低周波成分といった特定のものとわかっていれば、フィルタ処理でノイズを軽減させた波形に対してピーク検出を使うといった方法が有効です。

ローパスフィルタ、ハイパスフィルタ、バンドパスフィルタ、バンドストップフィルタについては、「信号処理」カテゴリのフィルタの項目を是非確認してみて下さい。

まとめ

本ページでは周波数波形におけるピークの重要性とピーク検出の概要、ノイズ処理の必要性を説明しました。

また、実際にPythonを使ったコードでサンプル波形のピーク検出を行い、正のピーク、負のピークが抽出されることと、ノイズの軽減に関するパラメータの有効性を確認しました。

実際の現場で使うためには、対象の波形の特徴を掴んだノイズ対策を行う必要がありますが、scipyを使えば簡単にピーク検出ができることがわかりました。

ピーク検出も信号処理の基本!またまた新しい武器を獲得しました!
Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています!

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

SNSでもご購読できます。

コメント

コメントを残す

*