Поиск пиков в данных

Использовать findpeaks для поиска значений и местоположений локальных максимумов в наборе данных.

Файл spots_num.mat содержит среднее количество солнечных пятен, наблюдаемых каждый год с 1749 по 2012 год. Данные доступны в НАСА.

Найдите максимумы и годы их появления. Постройте их график вместе с данными.

load('spots_num.mat')

[pks,locs] = findpeaks(avSpots);

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

Некоторые вершины очень близки друг к другу. Те, которые не повторяются через регулярные промежутки времени. За 50-летний период таких пиков примерно пять.

Чтобы лучше оценить продолжительность цикла, используйте findpeaks опять же, но на этот раз ограничить разделение пик-пик по меньшей мере шестью годами. Вычислите средний интервал между максимумами.

[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')

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)])

Преобразование Фурье действительно достигает пика на ожидаемой частоте, подтверждая 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')

Эти две оценки вполне совпадают.

См. также

dlmread | findpeaks

Связанные темы

• Пиковый анализ