findpeaks
局所的最大値
構文
説明
[___] = findpeaks(___, は、前の構文の入力引数に加えて、名前と値の引数を使用してオプションを指定します。Name=Value)
出力引数を指定しない findpeaks(___) は信号をプロットし、ピーク値を重ね合わせます。
例
3 つのピークをもつベクトルを定義します。
data = [25 8 15 5 6 10 10 3 1 20 7];
局所的最大値を求めます。findpeaks 関数は、ピークを出現順に返します。関数はピーク検出において 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; y = sum(Hgt.*(exp(-((x-Pos)./Wdt).^2)),1);
既定の設定で findpeaks を使用して信号のピークとそれらの位置を求めます。
[pks,locs] = findpeaks(y,x);
findpeaks を使用してピークをプロットします。各ピークにラベルを付けます。
findpeaks(y,x) text(locs+.02,pks,num2str((1:numel(pks))'))
最も高いピークから低いピークに並べ替えます。
[psor,lsor] = findpeaks(y,x,SortStr="descend");
findpeaks(y,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; y = sum(Hgt.*(exp(-((x-Pos)./Wdt).^2)),1) + base;
findpeaks を使用して 4 以上のプロミネンスをもつピークの場所を特定しプロットします。最も高いピークと最も低いピークのみが、最小プロミネンスの条件を満たします。
[~,locs4,~,proms4] = findpeaks(y,x,MinPeakProminence=4)
locs4 = 1×2
0.1982 0.4995
proms4 = 1×2
5.5773 4.4171
findpeaks(y,x,MinPeakProminence=4,Annotate="extents")
すべてのピークの位置、プロミネンス、およびプロミネンスの半分におけるピーク幅を表示します。2 番目および 4 番目のピークが、4 以上のプロミネンスをもつピークです。
[~,locs,widths,proms] = findpeaks(y,x)
locs = 1×6
0.1001 0.1982 0.2983 0.4995 0.7017 0.8018
widths = 1×6
0.0154 0.0431 0.0377 0.0625 0.0274 0.0409
proms = 1×6
2.6816 5.5773 3.1448 4.4171 2.9191 3.6363
黒点の平均個数データからピークを抽出し、平均黒点サイクルを推定します。
sunspot.dat ファイルを読み込みます。このファイルには、1700 年から 1987 年までに観測された太陽黒点数の毎年の平均が含まれています。年ごとの黒点数をプロットし、最大値を重ね合わせます。1855 年から 1920 年までの期間を拡大します。
load sunspot.dat year = sunspot(:,1); avSpots = sunspot(:,2); findpeaks(avSpots,year) xlabel("Year") xlim([1855 1920])
黒点は周期的現象です。その数は、おおよそ 11 年ごとにピークを迎えることが知られています。
ピーク間の平均間隔を求めることで、サイクルの期間を推定します。まず、最大値ではないピークを無視します。受け入れられるピーク間の間隔を、6 年を超える周期に制限します。
findpeaks(avSpots,year,MinPeakDistance=6)
xlabel("Year")
xlim([1855 1920])
年次データを使用して datetime 配列を作成します。黒点は毎年春分近くの 3 月 20 日にカウントされたと仮定します。黒点数がピークとなる年を求めます。関数 years を使用して、最小ピーク間隔を duration として指定します。
ty = datetime(year,3,20); [pk,lk] = findpeaks(avSpots,ty,MinPeakDistance=years(6)); plot(ty,avSpots,lk,pk,"o") xlabel("Year")
平均黒点サイクルを計算します。推定された平均黒点サイクルは、既知の 11 年周期とよく一致しています。
dttmCycle = years(mean(diff(lk)))
dttmCycle = 10.9600
特定の特徴条件 (最小または最大のピーク距離や高さ、プロミネンス、および隣接値との差のしきい値など) を満たすピークを求めることができます。findpeaks を使用して、1 つまたは複数のピーク特徴に制約をかけます。
音声信号の特徴によるピークの制約
7418 Hz でサンプリングされた音声信号を読み込みます。200 サンプルを選択します。
load mtlb
select = mtlb(1001:1200);音声信号からすべてのピークを検出し、強調表示します。次に、少なくとも 5 ms の間隔があるピークを検出し、強調表示します。最小間隔の制約を適用するために、findpeaks は信号内で最も高いピークを選択し、そのピークから 5 ms 以内にあるすべてのピークを除外します。次に、関数は残りのピークの中で最も高いピークを選択するためにこの手順を繰り返し、考慮すべきピークがなくなるまで反復します。
figure tiledlayout("vertical") nexttile findpeaks(select,Fs) legend("No Constraints",Location="southwest") nexttile lgn = @(txt) legend(txt,FontName="FixedWidth",Location="southwest"); findpeaks(select,Fs,MinPeakDistance=0.005) lgn("MinPeakDistance")
以下の各制約を個別に満たすピークを検出し、プロットします。
1 V 以上の振幅をもつピークを検出します。
隣接するサンプルよりも 1 V 以上高いピークを検出します。
より高い値に信号が達する前にいずれかの側に 1 V 以上ドロップしているピークを検出します。
figure tiledlayout("vertical") nexttile findpeaks(select,Fs,MinPeakHeight=1) lgn("MinPeakHeight") nexttile findpeaks(select,Fs,Threshold=1) lgn("Threshold") nexttile findpeaks(select,Fs,MinPeakProminence=1) lgn("MinPeakProminence")
チャープ信号による複数のピーク特徴の制約
1000 サンプルのチャープ信号から、幅が 2.5 ms から 4 ms の範囲にあるすべてのピークを検出し、強調表示します。
N = 1000; x = sin(2*pi*(1:N)/N + (10*(1:N)/N).^2); figure findpeaks(x,Fs,Annotate="extents", ... MinPeakWidth=2.5e-3,MaxPeakWidth=4e-3) title("Peaks with Multiple Constraints")
データが指定した飽和点より大きい場合、センサーはクリップした読み取り値を返します。これらの極値を無意味なものとして無視するか、あるいは解析に組み込むかを選択できます。
分散 0.1² のホワイト ガウス ノイズに含まれる、周波数が 5 Hz と 3 Hz の三角関数の積で構成された信号を生成します。100 Hz のレートで信号を 1 秒間サンプリングします。再現可能な結果が必要な場合は、乱数発生器をリセットします。
rng("default")
Fs = 100;
t = 0:1/Fs:1;
s = sin(2*pi*5*t).*sin(2*pi*3*t) + randn(size(t))/10;絶対値が 0.32 の指定範囲より大きいすべての読み取り値を切り捨てることで、飽和した測定値を合成します。次に、信号を正規化して 0.1 から 1 の範囲になるようにします。正規化した振幅をデシベルに変換し、結果をプロットします。
bnd = 0.32; s(s>bnd) = bnd; s(s<-bnd) = -bnd; sNorm = 0.1+0.9*(s-min(s))/(max(s)-min(s)); sig = mag2db(sNorm); plot(t,sig) xlabel("Time (s)") ylabel("Magnitude (dB)") ylim([-25 5])
信号のピークと谷を特定します。findpeaks は、平らなピークのそれぞれの立ち上がりエッジのみを識別します。
[pk,lcp] = findpeaks(sig,t); % Peaks [vl,lcv] = findpeaks(-sig,t); % Valleys hold on plot(lcp,pk,"x",lcv,-vl,"*")
名前と値の引数 Threshold を使用して、平らなピークと谷を除外します。ピークまたは谷とその隣接値との間には、最小振幅差として 0.1 dB が必要です。
[pkt,lctp] = findpeaks(sig,t,Threshold=0.1); [vlt,lctv] = findpeaks(-sig,t,Threshold=0.1); plot(lctp,pkt,"o",LineWidth=1.2) plot(lctv,-vlt,"d",LineWidth=1.2)
ガウス曲線の和を構成する信号を作成します。各曲線の位置、高さ、幅を指定します。
x = linspace(0,1,1000); Pos = [1 2 3 5 7 8]'/10; Hgt = [7 6 3 2 2 3]'; Wdt = [3 8 4 3 4 6]'/100; y = sum(Hgt.*(exp(-((x-Pos)./Wdt).^2)),1);
プロミネンスの半分および高さの半分を基準として使用し、ピークの幅を測定します。
tiledlayout("flow") nexttile findpeaks(y,x,Annotate="extents") title("Half-Prominence Peak Widths") nexttile findpeaks(y,x,Annotate="extents",WidthReference="halfheight") title("Half-Height Peak Widths")
"x" 軸で少なくとも 0.5 単位の間隔がある、最も高いピークを選択します。プロミネンスの半分および高さの半分を基準として使用し、ピークの幅を測定します。
figure tiledlayout("flow") nexttile findpeaks(y,x,MinPeakDistance=0.5,Annotate="extents") title("Half-Prominence Peak Widths") nexttile findpeaks(y,x,MinPeakDistance=0.5,Annotate="extents", ... WidthReference="halfheight") title("Half-Height Peak Widths")
最小間隔の条件を満たすのは最初と最後のピークのみであるため、プロットに表示されるピーク幅はこれら 2 つのピークに対応します。各ピークの範囲は変化しないため、指定された条件や、そのピークが選択されたかどうかにかかわらず、ピーク幅の値は保持されます。
sinc 関数カーネルを使った非線形最小二乗法を使用して、信号内の主要な 2 つのピークの正確なピーク位置推定値を取得します。
信号の生成
線形 FM 波形のレーダー パルス圧縮により、sinc 形状のスペクトルが生成されます。この際、ピークの周波数位置は、レーダーと検出されたオブジェクトとの間の距離に比例します。まず、findpeaksを使用してピークの位置と振幅を推定し、次に、refinepeaksを使用して推定結果を改善することができます。この例では、合成ノイズレス パルス圧縮信号のピーク振幅と位置を求め、refinepeaks を使用して推定結果を改善します。
それぞれ 4.76 kHz と 35.8 kHz のピークが 1 と 1.5 である 2 つの sinc 形状の波形で構成される信号を生成します。周波数間隔を 2.5 Hz に設定します。
aTg = [1 1.5]; fTg = 1e3*[4.76 35.8]; freqkHzFull = (0:0.0025:50)'; waveFull = abs(sinc([1 0.5].*(freqkHzFull-fTg/1e3)))*aTg';
信号を係数 200 でダウンサンプリングし、サンプル間の周波数間隔が 0.5 kHz になるようにします。この例では、ダウンサンプリングされた信号のピークの振幅と位置の推定値を改善し、改善後の推定値を元の信号の値と比較します。
freq = downsample(freqkHzFull,200); wave = downsample(waveFull,200); plot(freqkHzFull,waveFull,freq,wave,"*") legend(["Full signal" "Selected samples"],Location="northwest") xlabel("Frequency (kHz)") ylabel("Magnitude")
非線形最小二乗法を使用したピークの改善
findpeaks を使用して、信号に含まれる最も高い 2 つのピークの振幅、位置、および半値幅の初期推定値を求めます。
[PV,PL,PW] = findpeaks(wave,NPeaks=2, ... SortStr="descend",WidthReference="halfheight"); table(PV,freq(PL), RowNames="Peak estimate "+(1:numel(PV)), ... VariableNames=["Amplitude" "Frequency (kHz)"])
ans=2×2 table
Amplitude Frequency (kHz)
_________ _______________
Peak estimate 1 1.4824 36
Peak estimate 2 0.9374 5
refinepeaks を使用して、非線形最小二乗法 (NLS) でピーク推定結果を改善します。信号の周波数点とピーク幅を指定します。ピーク値は予想値の 1.5 と 1 にかなり近く、周波数位置はそれぞれ 35.8 kHz と 4.76 kHz によく近似しています。
LW = max(PW,2); [Ypk,Xpk] = refinepeaks(wave,PL,freq,Method="NLS",LobeWidth=LW); table(Ypk,Xpk, RowNames="Refined peak "+(1:numel(Ypk)), ... VariableNames=["Amplitude" "Frequency (kHz)"])
ans=2×2 table
Amplitude Frequency (kHz)
_________ _______________
Refined peak 1 1.5063 35.8
Refined peak 2 1.0163 4.7628
改善後のピークの振幅を y 軸に、初期ピーク推定値を基準とした更新後のピークの位置を x 軸にプロットします。最初に推定された 2 つのピークとその周囲の 2 つのサンプルは、それぞれ 0.5 kHz ずつ離れています。塗りつぶされた円で示される改善後のピークは、最初に推定されたピーク位置に対する実際のピーク位置と、修正された振幅を示しています。
refinepeaks(wave,PL,freq,Method="NLS",LobeWidth=LW) yline(aTg) % Theoretical peak amplitudes errorBounds = aTg.*(1+0.03*[-1;1]); yline(errorBounds(:),":") % ±3% error bounds legend("Peak "+[1 2])
入力引数
名前と値の引数
出力引数
詳細
ヒント
最初に findpeaks を使用して信号のピークを推定し、次に refinepeaks を使用してその振幅と位置を改善することができます。
振幅が y で位置が x である信号があると仮定します。次のコードの抜粋は、y と x からピークを推定して改善する方法を示しています。
[yPeaks,xPeaksIdx] = findpeaks(y); [yRPeaks,xRPeaks] = refinepeaks(y,xPeaksIdx,x)
