findpeaks

pks = findpeaks(data) は、入力信号ベクトル data から局所的最大値 (ピーク) をもつベクトルを返します。"局所的なピーク" は、2 つの隣接するサンプルよりも大きいか、Inf と等しいデータサンプルです。Inf でない信号の端点は除外されます。ピークが平坦である場合、関数は最小のインデックスをもつ点のみを返します。

[pks,locs] = findpeaks(data) はさらに、ピークが発生するインデックスを返します。

[pks,locs,w,p]= findpeaks(data) は、さらにベクトル w としてピーク幅を、ベクトル p としてピークのプロミネンスを返します。

[___] = findpeaks(data,x) は x を位置ベクトルとして指定し、前の構文の出力引数のいずれかを返します。locs および w はx で表されます。

[___] = findpeaks(data,Fs) は、データのサンプルレート Fs を指定します。data の最初のサンプルは時間ゼロでとられたものと想定されます。locs および w は時間の単位に変換されます。

[___] = findpeaks(___,Name,Value) は、前の構文の入力引数のいずれかと、名前と値のペアを使用してオプションを指定します。

出力引数を指定しない findpeaks(___) は信号をプロットし、ピーク値を重ね合わせます。

例

ベクトルにおけるピークを検出します。

3 つのピークをもつベクトルを定義し、プロットします。

data = [25 8 15 5 6 10 10 3 1 20 7];
plot(data)

局所的最大値を求めます。ピークは発生順に出力されます。最初のサンプルは最大であっても含まれません。平坦なピークの場合、関数は最も低いインデックスをもつ点のみを返します。

pks = findpeaks(data)

pks = 1×3

    15    10    20

出力引数なしで findpeaks を使用してピークを表示します。

findpeaks(data)

ピークとその位置の検索

ベル型曲線の和を構成する信号を作成します。各曲線の位置、高さ、幅を指定します。

x = linspace(0,1,1000);

Pos = [1 2 3 5 7 8]/10;
Hgt = [3 4 4 2 2 3];
Wdt = [2 6 3 3 4 6]/100;

for n = 1:length(Pos)
    Gauss(n,:) = Hgt(n)*exp(-((x - Pos(n))/Wdt(n)).^2);
end

PeakSig = sum(Gauss);

個々の曲線とその和をプロットします。

plot(x,Gauss,'--',x,PeakSig)

既定の設定で findpeaks を使用して信号のピークとそれらの位置を求めます。

[pks,locs] = findpeaks(PeakSig,x);

findpeaks を使用してピークをプロットし、それにラベルを付けます。

findpeaks(PeakSig,x)

text(locs+.02,pks,num2str((1:numel(pks))'))

最も高いピークから低いピークに並べ替えます。

[psor,lsor] = findpeaks(PeakSig,x,'SortStr','descend');

findpeaks(PeakSig,x)

text(lsor+.02,psor,num2str((1:numel(psor))'))

ピークのプロミネンス

余弦 1 周期の区間を進むベル型曲線の和から構成される信号を作成します。各曲線の位置、高さ、幅を指定します。

x = linspace(0,1,1000);

base = 4*cos(2*pi*x);

Pos = [1 2 3 5 7 8]/10;
Hgt = [3 7 5 5 4 5];
Wdt = [1 3 3 4 2 3]/100;

for n = 1:length(Pos)
    Gauss(n,:) =  Hgt(n)*exp(-((x - Pos(n))/Wdt(n)).^2);
end

PeakSig = sum(Gauss)+base;

個々の曲線とその和をプロットします。

plot(x,Gauss,'--',x,PeakSig,x,base)

findpeaks を使用して 4 以上のプロミネンスをもつピークの場所を特定しプロットします。

findpeaks(PeakSig,x,'MinPeakProminence',4,'Annotate','extents')

最も高いピークと最も低いピークだけが条件を満たすピークです。

すべてのピークのプロミネンスおよびプロミネンスの半分の幅を表示します。

[pks,locs,widths,proms] = findpeaks(PeakSig,x);
widths

widths = 1×6

    0.0154    0.0431    0.0377    0.0625    0.0274    0.0409

proms

proms = 1×6

    2.6816    5.5773    3.1448    4.4171    2.9191    3.6363

最小離隔距離のピークの検出

黒点は周期的現象です。その数は、おおよそ 11 年ごとにピークを迎えることが知られています。

sunspot.dat ファイルを読み込みます。このファイルには、1700 年から 1987 年までに観測された太陽黒点数の毎年の平均が含まれています。最大値を求めてプロットします。

load sunspot.dat

year = sunspot(:,1);
avSpots = sunspot(:,2);

findpeaks(avSpots,year)

周期期間の推定は、相互の距離が非常に近いピークを無視することで精度が上がります。再度ピークを求めてプロットしますが、今回は受け入れることができるピーク間の離隔距離を 6 年よりも大きな値に制限します。

findpeaks(avSpots,year,'MinPeakDistance',6)

findpeaks で返されたピークの位置を使用して最大値間の平均間隔を計算します。

[pks,locs] = findpeaks(avSpots,year,'MinPeakDistance',6);

meanCycle = mean(diff(locs))

meanCycle = 10.9600

年データを使用して datetime 配列を作成します。黒点は毎年春分近くの 3 月 20 日にカウントされたと仮定します。黒点数がピークとなる年を求めます。関数 years を使用して、最小ピーク離隔距離を duration として指定します。

ty = datetime(year,3,20);

[pk,lk] = findpeaks(avSpots,ty,'MinPeakDistance',years(6));

plot(ty,avSpots,lk,pk,'o')

datetime の機能を使用して、平均黒点周期を計算します。

dttmCycle = years(mean(diff(lk)))

dttmCycle = 10.9600

データを使用して timetable を作成します。時間変数を年で指定します。データをプロットします。timetable の最後の 5 つのエントリを表示します。

TT = timetable(years(year),avSpots);
plot(TT.Time,TT.Variables)

entries = TT(end-4:end,:)

entries=5×1 timetable
      Time      avSpots
    ________    _______

    1983 yrs     66.6  
    1984 yrs     45.9  
    1985 yrs     17.9  
    1986 yrs     13.4  
    1987 yrs     29.3

ピーク特性の制約

7418 Hz でサンプリングされた音声信号を読み込みます。200 サンプルを選択します。

load mtlb
select = mtlb(1001:1200);

最小分離間隔 5 ms でピークを検出します。

この制約を適用するため、findpeaks は信号の最も高いピークを選択し、5 ms 以内のすべてのピークを排除します。次に、関数は残りのピークの中で最も高いピークを選択するためにこの手順を繰り返し、考慮すべきピークがなくなるまで反復します。

findpeaks(select,Fs,'MinPeakDistance',0.005)

1 V 以上の振幅をもつピークを検出します。

findpeaks(select,Fs,'MinPeakHeight',1)

隣接するサンプルよりも 1 V 以上高いピークを検出します。

findpeaks(select,Fs,'Threshold',1)

より高い値に信号が達する前にいずれかの側に 1 V 以上ドロップしているピークを検出します。

findpeaks(select,Fs,'MinPeakProminence',1)

飽和した信号のピーク

データが指定した飽和点より大きい場合、センサーはクリップした読み取り値を返します。これらのピークを無意味なものとして無視するか、あるいは解析に組み込むように選択できます。

分散 0.1² のホワイトガウスノイズに含まれる、周波数が 5 Hz と 3 Hz の三角関数の積で構成された信号を生成します。信号は、100 Hz のレートで 1 秒間サンプリングされます。再現可能な結果が必要な場合は、乱数発生器をリセットします。

rng default

fs = 1e2;
t = 0:1/fs:1-1/fs;

s = sin(2*pi*5*t).*sin(2*pi*3*t)+randn(size(t))/10;

0.32 の指定範囲より大きいすべての読み取り値を切り捨てることで、飽和した測定値をシミュレートします。飽和した信号をプロットします。

bnd = 0.32;
s(s>bnd) = bnd;

plot(t,s)
xlabel('Time (s)')

信号のピークの位置を特定します。findpeaks は、平らなピークのそれぞれの立ち上がりエッジのみをレポートします。

[pk,lc] = findpeaks(s,t);

hold on
plot(lc,pk,'x')

'Threshold' の名前と値のペアを使用して平らなピークを除外します。ピークとその両隣との間には、最小振幅差として $1 0^{- 4}$ が必要です。

[pkt,lct] = findpeaks(s,t,'Threshold',1e-4);

plot(lct,pkt,'o','MarkerSize',12)

ピーク幅の特定

ベル型曲線の和を構成する信号を作成します。各曲線の位置、高さ、幅を指定します。

x = linspace(0,1,1000);

Pos = [1 2 3 5 7 8]/10;
Hgt = [4 4 2 2 2 3];
Wdt = [3 8 4 3 4 6]/100;

for n = 1:length(Pos)
    Gauss(n,:) =  Hgt(n)*exp(-((x - Pos(n))/Wdt(n)).^2);
end

PeakSig = sum(Gauss);

個々の曲線とその和をプロットします。

plot(x,Gauss,'--',x,PeakSig)
grid

プロミネンスの半分を基準として使用してピーク幅を測定します。

findpeaks(PeakSig,x,'Annotate','extents')

今度は基準として高さの半分を使用して幅を再測定します。

findpeaks(PeakSig,x,'Annotate','extents','WidthReference','halfheight')
title('Signal Peak Widths')

入力引数

`data` — 入力データ
ベクトル

ベクトルとして指定される入力データ。data は実数で要素数が 3 以上でなければなりません。

データ型: double | single

`x` — 位置
ベクトル | `datetime` 配列

位置。ベクトルまたは datetime 配列で指定されます。x は単調増加で data と同じ長さでなければなりません。x を省略した場合、data のインデックスが位置として使用されます。

データ型: double | single | datetime

`Fs` — サンプルレート
正のスカラー

サンプルレート。正のスカラーで指定します。サンプルレートは単位時間あたりのサンプル数です。時間の単位が秒の場合、サンプルレートの単位は Hz です。

データ型: double | single

名前と値のペアの引数

オプションの引数 Name,Value のコンマ区切りペアを指定します。Name は引数名で、Value は対応する値です。Name は引用符で囲まなければなりません。Name1,Value1,...,NameN,ValueN のように、複数の名前と値のペアの引数を、任意の順番で指定できます。

例: 'SortStr','descend','NPeaks',3 では、信号の最も高いピークを 3 つ求めます。

`'NPeaks'` — ピークの最大数
正の整数スカラー

返されるピークの最大数。'NPeaks' および正の整数スカラーで構成されるコンマ区切りのペアとして指定します。findpeaks は入力データの最初の要素から演算を開始し、ピークの数が 'NPeaks' の値に達すると終了します。

データ型: double | single

`'SortStr'` — ピークの並べ替え
`'none'` (既定値) | `'ascend'` | `'descend'`

ピークの並べ替え。コンマ区切りの 'SortStr' と次の値のいずれかとのペアとして指定します。

'none' は入力データの発生順でピークを返します。
'ascend' はピークを最小値から最大値への昇順、つまり増加順で返します。
'descend' はピークを最大値から最小値への降順で返します。

`'MinPeakHeight'` — ピーク高の最小値
`-Inf` (既定値) | 実数スカラー

ピーク高の最小値。'MinPeakHeight' および実数スカラーで構成されるコンマ区切りペアとして指定します。この引数を使用して、findpeaks が 'MinPeakHeight' よりも高いピークのみを返すようにします。最小ピーク高を指定することで処理時間が短縮されます。

データ型: double | single

`'MinPeakProminence'` — 最小ピークのプロミネンス
0 (既定値) | 実数スカラー

最小ピークのプロミネンス。'MinPeakProminence' および実数スカラーで構成されるコンマ区切りペアとして指定します。この引数を使用して、findpeaks が少なくとも 'MinPeakProminence' の相対的な重要性をもつピークのみを返すようにします。詳細は、プロミネンスを参照してください。

データ型: double | single

`'Threshold'` — 最小高さ距離
0 (既定値) | 非負の実数スカラー

ピークとその両隣接値との高さの最小差異。'Threshold' および非負の実数スカラーで構成されるコンマ区切りペアとして指定します。findpeaks はこの引数を使用して、隣接値を少なくとも 'Threshold' の値分を上回るピークのみを返します。

データ型: double | single

`'MinPeakDistance'` — 最小ピーク離隔距離
0 (既定値) | 正の実数のスカラー

最小ピーク離隔距離。'MinPeakDistance' と正の実数のスカラーで構成されるコンマ区切りペアとして指定します。'MinPeakDistance' に値を指定した場合、アルゴリズムは信号の最も高いピークを選択し、信号の 'MinPeakDistance' 以内のすべてのピークを無視します。次に、関数は残りのピークの中で最も高いピークを選択するためにこの手順を繰り返し、考慮すべきピークがなくなるまで反復します。

位置ベクトル x を指定した場合、'MinPeakDistance' は必ず x で表されます。x が datetime 配列の場合は、'MinPeakDistance' を duration スカラーまたは日単位で表される数値スカラーとして指定します。
サンプルレート Fs を指定した場合、'MinPeakDistance' は必ず時間の単位で表されます。
x と Fs のどちらも指定しない場合、'MinPeakDistance' は必ずサンプルの単位で表されます。

この引数を使用して、findpeaks により大きなピークの近傍で発生する小さなピークを無視させます。

データ型: double | single | duration

`'WidthReference'` — 幅の測定値の基準の高さ
`'halfprom'` (既定値) | `'halfheight'`

幅の測定値の基準の高さ。'WidthReference' および 'halfprom' または 'halfheight' のいずれかで構成されるコンマ区切りペアとして指定します。findpeaks は降順の信号が水平の基準線を横切る点の間の距離としてピーク幅を推定します。線の高さは 'WidthReference' で指定されている基準を使用して選択されます。

'halfprom' はピークのプロミネンスの半分に等しい垂直距離にあるピークの下に基準線を配置します。詳細は、プロミネンスを参照してください。
'halfheight' はピーク高さの 1/2 の位置に基準線を配置します。この線は、任意のインターセプトポイントが 'MinPeakHeight'、'MinPeakProminence' および 'Threshold' の設定によって選択されたピークの境界を超えた場合に打ち切られます。ピーク間の境界は、ピーク間の最も低い谷の水平位置で定義されます。高さがゼロ未満のピークは破棄されます。

インターセプトポイントの位置は、線形内挿で計算されます。

`'MinPeakWidth'` — 最小ピーク幅
0 (既定値) | 正の実数のスカラー

最小ピーク幅。'MinPeakWidth' と正の実数のスカラーで構成されるコンマ区切りペアとして指定します。この引数を使用して、'MinPeakWidth' 以上の幅をもつピークだけを選択します。

位置ベクトル x を指定した場合、'MinPeakWidth' は必ず x で表されます。x が datetime 配列の場合は、'MinPeakWidth' を duration スカラーまたは日単位で表される数値スカラーとして指定します。
サンプルレート Fs を指定した場合、'MinPeakWidth' は必ず時間の単位で表されます。
x と Fs のどちらも指定しない場合、'MinPeakWidth' は必ずサンプルの単位で表されます。

データ型: double | single | duration

`'MaxPeakWidth'` — 最大ピーク幅
`Inf` (既定値) | 正の実数のスカラー

最大ピーク幅。'MaxPeakWidth' と正の実数のスカラーで構成されるコンマ区切りペアとして指定します。この引数を使用して、最大 'MaxPeakWidth' の幅をもつピークだけを選択します。

位置ベクトル x を指定した場合、'MaxPeakWidth' は必ず x で表されます。x が datetime 配列の場合は、'MaxPeakWidth' を duration スカラーまたは日単位で表される数値スカラーとして指定します。
サンプルレート Fs を指定した場合、'MaxPeakWidth' は必ず時間の単位で表されます。
x と Fs のどちらも指定しない場合、'MaxPeakWidth' は必ずサンプルの単位で表されます。

データ型: double | single | duration

`'Annotate'` — プロットのスタイル
`'peaks'` (既定値) | `'extents'`

プロットのスタイル。'Annotate' と次の値のいずれかで構成されるコンマ区切りペアとして指定します。

'peaks' は信号をプロットし、各ピークの位置と値の注釈を付けます。
'extents' は信号をプロットし、各ピークの位置、値、幅およびプロミネンスに注釈を付けます。

この引数は出力引数で findpeaks を呼び出している場合に無視されます。

出力引数

`pks` — 局所的最大値
ベクトル

局所的最大値。信号値のベクトルで返されます。局所的最大値がない場合、pks は空です。

`locs` — ピークの位置
ベクトル

ベクトルとして返されるピークの位置。

位置ベクトル x を指定した場合、locs はピークのインデックスに x の値を含めます。
サンプルレート Fs を指定した場合、locs は連続するサンプル間の時間差が 1/Fs となる時点の数値ベクトルです。
x と Fs のどちらも指定しない場合、locs は整数インデックスのベクトルです。

`w` — ピーク幅
ベクトル

ピーク幅。実数のベクトルとして返します。各ピークの幅は、WidthReference で指定された高さの基準線に信号が交差する、ピークの左右の点の間の距離を計算します。点自体は線形内挿で検出されます。

位置ベクトル x を指定した場合、幅は x で表されます。
サンプルレート Fs を指定した場合、時間の単位で表されます。
x と Fs のどちらも指定しない場合、幅はサンプル単位で表されます。

`p` — ピークのプロミネンス
ベクトル

ピークのプロミネンス。実数のベクトルとして返します。ピークのプロミネンスとは、信号がピークよりも高いレベルに戻るか、端点に到達する前にピークのいずれかの側に信号が下降して到達する必要のある最小垂直距離です。詳細は、プロミネンスを参照してください。

詳細