Использование 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
Некоторый peaks очень близки друг к другу. Те таковые, которые не повторяются при регулярных интервалах. На 50-летний период такого peaks примерно пять.
Чтобы лучше оценить длительность цикла, используйте 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')
Две оценки совпадают довольно хорошо.