Найдите Peaks в данных
Используйте findpeaks, чтобы найти значения и местоположения локальных максимумов в наборе данных.
Файл spots_num.mat содержит среднее количество солнечных пятен, наблюдаемых каждый год от 1 749 до 2012. Данные доступны от НАСА.
Найдите максимумы и их годы вхождения. Постройте их наряду с данными.
load(fullfile(matlabroot,'examples','signal','spots_num.mat')) [pks,locs] = findpeaks(avSpots); plot(year,avSpots,year(locs),pks,'or') xlabel('Year') ylabel('Number') axis tight
Некоторый peaks очень друг близко к другу. Те, которые не являются, повторяются равномерно. Существует примерно пять таких peaks на 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')
Две оценки совпадают вполне хорошо.