Анализ циклических данных с помощью БПФ

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

За почти 300 лет астрономы свели в таблицу количество и размер солнечных пятен с помощью относительного числа солнечных пятен Цюриха. Постройте график цюрихского числа приблизительно за 1700 - 2000 годы.

load sunspot.dat
year = sunspot(:,1);
relNums = sunspot(:,2);
plot(year,relNums)
xlabel('Year')
ylabel('Zurich Number')
title('Sunspot Data')

Figure contains an axes. The axes with title Sunspot Data contains an object of type line.

Чтобы присмотреться к циклической природе активности солнечных пятен, постройте график первых 50 лет данных.

plot(year(1:50),relNums(1:50),'b.-');
xlabel('Year')
ylabel('Zurich Number')
title('Sunspot Data')

Figure contains an axes. The axes with title Sunspot Data contains an object of type line.

Преобразование Фурье является фундаментальным инструментом в обработке сигналов, который идентифицирует частотные составляющие в данных. Использование fft function, взять преобразование Фурье данных Цюриха. Удалите первый элемент выхода, который хранит сумму данных. Постройте график оставшейся части выхода, который содержит зеркальное изображение комплексных коэффициентов Фурье вокруг действительной оси.

y = fft(relNums);
y(1) = [];
plot(y,'ro')
xlabel('real(y)')
ylabel('imag(y)')
title('Fourier Coefficients')

Figure contains an axes. The axes with title Fourier Coefficients contains an object of type line.

Коэффициенты Фурье сами по себе трудно интерпретировать. Более значимой мерой коэффициентов является их величина в квадрате, которая является мерой степени. Поскольку половина коэффициентов повторяется в величине, вам нужно только вычислить степень на одной половине коэффициентов. Постройте график спектра степени как функции от частоты, измеряемой в циклах в год.

n = length(y);
power = abs(y(1:floor(n/2))).^2; % power of first half of transform data
maxfreq = 1/2;                   % maximum frequency
freq = (1:n/2)/(n/2)*maxfreq;    % equally spaced frequency grid
plot(freq,power)
xlabel('Cycles/Year')
ylabel('Power')

Figure contains an axes. The axes contains an object of type line.

Максимальная активность солнечных пятен происходит реже, чем один раз в год. Для представления о циклической активности, которую легче интерпретировать, постройте график степени как функции от периода, измеряемой в годах за цикл. График показывает, что активность солнечных пятен достигает пика примерно раз в 11 лет.

period = 1./freq;
plot(period,power);
xlim([0 50]); %zoom in on max power
xlabel('Years/Cycle')
ylabel('Power')

Figure contains an axes. The axes contains an object of type line.

См. также

| |

Похожие темы