Python/SciPy のピーク検出アルゴリズム
-
19-09-2019 - |
質問
一次導関数のゼロクロッシングなどを見つけて自分で何かを書くこともできますが、標準ライブラリに含めるのに十分な一般的な関数のように思えます。誰か知っていますか?
私の特定のアプリケーションは 2D 配列ですが、通常は FFT などでピークを見つけるために使用されます。
具体的には、この種の問題では、複数の強いピークがあり、その後、無視すべきノイズによって引き起こされる小さな「ピーク」が多数存在します。これらは単なる例です。私の実際のデータではありません:
1次元ピーク:
2次元ピーク:
ピーク検出アルゴリズムは、これらのピーク (値だけでなく) の位置を見つけ、理想的には、最大値を持つインデックスだけでなく、真のサンプル間のピークを見つけるでしょう。 二次補間 か何か。
通常、関心があるのは少数の強いピークだけであるため、それらは特定のしきい値を超えているか、最初のピークであるため選択されます。 n 振幅によってランク付けされた、順序付きリストのピーク。
先ほども言いましたが、私はこのようなものを自分で書く方法を知っています。私は、うまく動作することが知られている既存の関数またはパッケージがあるかどうかを尋ねているだけです。
アップデート:
私 MATLAB スクリプトを翻訳しました 1-D の場合は適切に機能しますが、さらに良くなる可能性があります。
更新されたアップデート:
シックステンベ より良いバージョンを作成しました 1-D の場合。
解決
あなたが探しているものはSciPyによって提供されているとは思えません。この状況では、私なら自分でコードを書きます。
scipy.interpolate によるスプライン補間と平滑化は非常に優れており、ピークをフィッティングして最大値の位置を見つけるのに非常に役立ちます。
他のヒント
私は同様の問題を検討していますが、最良の参考文献のいくつかは化学からのものであることがわかりました(質量分析データでのピーク検出から)。ピーキング検出アルゴリズムの詳細なレビューについては、以下をお読みください。 これ. 。これは、私がこれまでに出会ったピーク発見テクニックに関する最も明確なレビューの 1 つです。(ウェーブレットは、ノイズの多いデータでこの種のピークを見つけるのに最適です。)
ピークが明確に定義されており、ノイズに隠れていないように見えます。そのため、スムーズな savtizky-golay 導関数を使用してピークを見つけることをお勧めします (上記のデータを微分するだけでは、誤検知が大量に発生することになります)。これは非常に効果的な手法であり、実装も非常に簡単です (基本的な演算を備えた行列クラスが必要です)。最初の S-G 導関数のゼロクロスを見つけるだけで満足できると思います。
関数 scipy.signal.find_peaks
, は、その名前が示すように、これに役立ちます。ただし、パラメータをよく理解することが重要です width
, threshold
, distance
そして何よりも prominence
良好なピーク抽出を実現します。
私のテストとドキュメントによると、次の概念は プロミネンス は、良好なピークを維持し、ノイズの多いピークを破棄するための「便利な概念」です。
とは (地形的な)隆起?それは 「頂上からより高い地形に行くために下降するのに必要な最低高さ」, ここでわかるように、
アイデアは次のとおりです。
プロミネンスが高いほど、ピークはより「重要」になります。
テスト:
多くの困難を伴うため、(ノイズの多い) 周波数が変化する正弦波を意図的に使用しました。ということがわかります。 width
最小値を設定した場合、パラメータはここではあまり役に立ちません。 width
高すぎると、高周波部分の非常に近いピークを追跡できなくなります。設定した場合 width
低すぎると、信号の左側に不要なピークが多数発生します。同じ問題 distance
. threshold
直接の近隣とのみ比較しますが、ここでは役に立ちません。 prominence
最良の解決策を提供するものです。これらのパラメータの多くは組み合わせられることに注意してください。
コード:
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import find_peaks
x = np.sin(2*np.pi*(2**np.linspace(2,10,1000))*np.arange(1000)/48000) + np.random.normal(0, 1, 1000) * 0.15
peaks, _ = find_peaks(x, distance=20)
peaks2, _ = find_peaks(x, prominence=1) # BEST!
peaks3, _ = find_peaks(x, width=20)
peaks4, _ = find_peaks(x, threshold=0.4) # Required vertical distance to its direct neighbouring samples, pretty useless
plt.subplot(2, 2, 1)
plt.plot(peaks, x[peaks], "xr"); plt.plot(x); plt.legend(['distance'])
plt.subplot(2, 2, 2)
plt.plot(peaks2, x[peaks2], "ob"); plt.plot(x); plt.legend(['prominence'])
plt.subplot(2, 2, 3)
plt.plot(peaks3, x[peaks3], "vg"); plt.plot(x); plt.legend(['width'])
plt.subplot(2, 2, 4)
plt.plot(peaks4, x[peaks4], "xk"); plt.plot(x); plt.legend(['threshold'])
plt.show()
scipyにはという名前の関数があります scipy.signal.find_peaks_cwt
これはあなたのニーズに適しているように思えますが、私には経験がないので、お勧めできません。
http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.find_peaks_cwt.html
Python でどのピーク検出アルゴリズムを使用すればよいかわからない人のために、ここでは代替アルゴリズムの簡単な概要を示します。 https://github.com/MonsieurV/py-findpeaks
自分自身を MatLab と同等にしたい findpeaks
関数、私はそれを発見しました ピーク検出関数 マルコス・ドゥアルテのは良いキャッチです。
使い方はとても簡単です:
import numpy as np
from vector import vector, plot_peaks
from libs import detect_peaks
print('Detect peaks with minimum height and distance filters.')
indexes = detect_peaks.detect_peaks(vector, mph=7, mpd=2)
print('Peaks are: %s' % (indexes))
これにより、次のことが得られます。
信頼性の高い方法でスペクトル内のピークを検出することは、たとえば 80 年代の音楽/オーディオ信号の正弦波モデリングに関する研究など、かなり研究されてきました。文献で「正弦波モデリング」を探してください。
信号が例のようにクリーンであれば、単純に「N 個の近傍よりも振幅が大きいものをください」で十分にうまく機能するはずです。ノイズの多い信号がある場合、シンプルですが効果的な方法は、時間の経過とともにピークを確認し、追跡することです。次に、スペクトル ピークではなくスペクトル線を検出します。そうですね、信号のスライディング ウィンドウで FFT を計算して、時間内のスペクトルのセット (スペクトログラムとも呼ばれます) を取得します。次に、時間の経過に伴うスペクトル ピークの変化を調べます (つまり、連続したウィンドウで)。
データの外れ値を見つけるための標準的な統計関数とメソッドがあり、おそらく最初のケースで必要なものです。導関数を使用すると、2 つ目の問題が解決されます。ただし、連続関数とサンプリングデータの両方を解決する方法はわかりません。
まず最初に、さらなる仕様がなければ、「ピーク」の定義は曖昧です。たとえば、次のシリーズの場合、5-4-5 を 1 つのピークと呼びますか、それとも 2 つのピークと呼びますか?
1-2-1-2-1-1-5-4-5-1-1-5-1
この場合、少なくとも 2 つのしきい値が必要になります。1) 高いしきい値。このしきい値を超えると、極値がピークとして記録されます。2) 低いしきい値なので、それ以下の小さな値で区切られた極値が 2 つのピークになります。
ピーク検出は、極値理論の文献でよく研究されているトピックであり、「極値のデクラスタリング」としても知られています。その典型的な用途には、環境変数の継続的な読み取りに基づいて危険事象を識別することが含まれます。風速を分析して嵐のイベントを検出します。