Python でピークを見つける

局所最大値とは、近傍のサンプルより大きい値です。系列全体で最も大きい値とは限りません。1次元系列では、各サンプルを隣接する値と比較すればシンプルに検出できます。候補を条件で絞り込みたい場合は、scipy.signal.find_peaks()を使えます。これらの関数が返すのは元データに対する0始まりのインデックスなので、そのインデックスを使ってピーク値を取得します。

NumPyでシンプルな局所最大値を検出する

次の例では1次元配列を使い、検出した各ピークのインデックスと値を表示します。

import numpy as np

data = np.array([1, 3, 7, 1, 2, 6, 0, 5, 0, 4])
local_maxima = (data[1:-1] > data[:-2]) & (data[1:-1] > data[2:])
indices = np.flatnonzero(local_maxima) + 1
print("indices:", indices.tolist())
print("values:", data[indices].tolist())

indices: [2, 5, 7]
values: [7, 6, 5]

2つの比較により、サンプルが左右の隣接値の両方を厳密に上回る場合だけ検出されます。スライスで比較するのは内部のサンプルだけなので、先頭と末尾は候補になりません。同じ値が連続する平坦な部分も除外されます。インデックスに加えるオフセットは、短くなった比較配列の位置をdata内の位置に対応させるためのものです。この厳密な条件がデータに合う場合は、NumPyを使う方法が適しています。この入力で検出されるインデックスは[2, 5, 7]で、対応する値は[7, 6, 5]です。

scipy.signal.find_peaks()でピークを条件指定して絞り込む

一般的な1次元ピーク検出には、SciPyのscipy.signal.find_peaks()を使えます。この関数は局所最大値のインデックスを返し、任意の条件を適用できます。次の例では高さが3以上のピークを選び、インデックスと値の両方を表示します。

import numpy as np
from scipy.signal import find_peaks

data = np.array([1, 3, 7, 1, 2, 6, 0, 5, 0, 4])
indices, properties = find_peaks(data, height=3, distance=2, prominence=1)
print("indices:", indices.tolist())
print("values:", data[indices].tolist())

indices: [2, 5, 7]
values: [7, 6, 5]

heightはピークの値で絞り込み、distanceは残すピーク同士の間隔をサンプル数で指定します。prominenceはピークが周囲からどの程度際立っているかを表します。distanceは採用されるピークに影響しますが、プロパティとして返される値ではありません。ほかに、widthやthresholdなどの任意条件もあります。SciPyはインデックスとプロパティの辞書を返します。この辞書には、指定した条件に応じて得られたプロパティが含まれます。この入力と引数の場合、返されるインデックスは[2, 5, 7]、対応する値は[7, 6, 5]です。すべての局所最大値に意味があると考えるのではなく、信号の性質に合った条件を選びましょう。

scipy.signal.argrelextrema()で相対的な極値を検出する

指定した近傍のサンプルと比較したい場合は、scipy.signal.argrelextrema()を使います。最大値を探すにはnp.greaterなどの比較関数を渡し、比較対象に含める左右それぞれのサンプル数をorderで指定します。戻り値は入力の各次元に対応するインデックス配列のタプルです。1次元配列の場合、インデックスは[0]で取り出します。

import numpy as np
from scipy.signal import argrelextrema

data = np.array([1, 3, 7, 1, 2, 6, 0, 5, 0, 4])
indices = argrelextrema(data, np.greater, order=1)[0]
print("indices:", indices.tolist())
print("values:", data[indices].tolist())

indices: [2, 5, 7]
values: [7, 6, 5]

order=1では、各候補を左右それぞれ1つのサンプルと比較します。orderを大きくすると比較する近傍が広がり、検出される相対極値が変わることがあります。この入力でorder=1を指定した場合、返されるインデックスは[2, 5, 7]で、対応する値は[7, 6, 5]です。

代替手段としてdetecta.detect_peaks()を使う

detecta.detect_peaks()は、1次元系列向けのサードパーティ製の選択肢です。プロジェクトの環境にpython -m pip install detectaで依存パッケージをインストールしてから、サンプルデータをdetect_peaks()に渡します。

data = [1, 3, 7, 1, 2, 6, 0, 5, 0, 4]
try:
    from detecta import detect_peaks
except ImportError:
    print("Install detecta with: python -m pip install detecta")
else:
    indices = detect_peaks(data)
    values = [data[index] for index in indices]
    print("indices:", list(indices))
    print("values:", values)

サンプル系列で期待されるピークの0始まりのインデックスは[2, 5, 7]で、対応する値は[7, 6, 5]です。提示された実行証拠で確認できるのは依存パッケージがない場合のメッセージだけなので、このdetectaの結果は実行で確認された出力ではなく、比較用の期待値です。インポート時の確認処理は残してください。detectaを利用できない場合、検出を試みずにインストールコマンドを表示します。

ピーク検出方法を選ぶ

隣接する値を厳密に上回ることだけが必要な条件で、処理の仕組みも直接確認したい場合は、NumPyでの比較が適しています。高さ、間隔、顕著度、幅などの条件による絞り込みが必要なら、scipy.signal.find_peaks()を選びましょう。極値の定義で比較する近傍や比較関数が重要なら、scipy.signal.argrelextrema()を使います。サードパーティ製APIがプロジェクトに合い、依存パッケージを追加できる場合は、detecta.detect_peaks()も候補になります。

どの方法も、入力系列に対する0始まりのインデックスを返します。結果を解釈する前に、各方法の定義とデータが合っているか確認してください。端点の扱い、平坦な部分の扱い、絞り込み条件は方法によって異なります。ここで扱った例は1次元の数値データの局所最大値を対象としています。多次元データのピーク検出には、そのデータ形状に対応した方法が必要です。