Main Content

このページの翻訳は最新ではありません。ここをクリックして、英語の最新版を参照してください。

データのピークの検出

findpeaks を使用して、データセットから局所的最大値とその位置を検出します。

ファイル spots_num.mat には、1749 年から 2012 年までに毎年観測された太陽黒点数の平均が含まれています。このデータは NASA から入手できます。

最大値とその発生年を検出します。それらとデータをプロットします。

load('spots_num.mat')

[pks,locs] = findpeaks(avSpots);

plot(year,avSpots,year(locs),pks,'or')
xlabel('Year')
ylabel('Number')
axis tight

Figure contains an axes. The axes contains 2 objects of type line.

距離が非常に近いピークがいくつかあります。それらは定期的には繰り返されない種類のものです。そのようなピークはおおよそ 50 年に 1 つ存在し、5 個ほどあります。

周期の期間をより正確に推定するには、もう一度 findpeaks を使用します。ただし、今回はピーク間の間隔を 6 年以上に制限します。最大値間の平均間隔を計算します。

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

plot(year,avSpots,year(locs),pks,'or')
xlabel('Year')
ylabel('Number')
title('Sunspots')
axis tight

legend('Data','peaks','Location','NorthWest')

Figure contains an axes. The axes with title Sunspots contains 2 objects of type line. These objects represent Data, peaks.

cycles = diff(locs);
meanCycle = mean(cycles)
meanCycle = 10.8696

太陽の活動が約 11 年ごとに繰り返されていることはよく知られています。フーリエ変換を使って確認します。信号の変動に集中できるよう、平均値を減算します。サンプル レートが年単位で測定されていることを思い出してください。ナイキスト周波数以下の周波数を使用します。

Fs = 1;
Nf = 512;

df = Fs/Nf;
f = 0:df:Fs/2-df;

trSpots = fftshift(fft(avSpots-mean(avSpots),Nf));

dBspots = 20*log10(abs(trSpots(Nf/2+1:Nf)));

yaxis = [20 85];
plot(f,dBspots,1./[meanCycle meanCycle],yaxis)
xlabel('Frequency (year^{-1})')
ylabel('| FFT | (dB)')
axis([0 1/2 yaxis])
text(1/meanCycle + .02,25,['<== 1/' num2str(meanCycle)])

Figure contains an axes. The axes contains 3 objects of type line, text.

フーリエ変換はまさに予想していた周波数でピークとなり、11 年ごとという予想が確認されました。フーリエ変換の最も高いピークを特定して周期を求めることもできます。

[pk,MaxFreq] = findpeaks(dBspots,'NPeaks',1,'SortStr','descend');

Period = 1/f(MaxFreq)
Period = 10.8936
hold on
plot(f(MaxFreq),pk,'or')
hold off
legend('Fourier transform','1/meanCycle','1/Period')

Figure contains an axes. The axes contains 4 objects of type line, text. These objects represent Fourier transform, 1/meanCycle, 1/Period.

この 2 つの推定はぴったり一致します。

参考

|

関連するトピック