Используйте 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
Некоторые вершины расположены очень близко друг к другу. Те, которые не повторяются через регулярные промежутки времени. На 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)
Хорошо известно, что циклы солнечной активности происходят примерно каждые 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)
hold on
plot(f(MaxFreq),pk,'or')
hold off
legend('Fourier transform','1/meanCycle','1/Period')
Эти две оценки вполне совпадают.