Изучаем цикличные данные при помощи БПФ (FFT)¶
В этой демонстрации мы проанализируем данные, содержащие хронику наблюдения за количеством пятен на поверхности Солнца, и выявим циклические составляющие при помощи алгоритма быстрого преобразования Фурье.
Введение¶
Преобразование Фурье позволяет нам анализировать характер изменения данных, и конечно раньше всего они применялись для анализа физических и природных процессов.
Данные, которые мы рассмотрим, являются сборником почти 300-летней истории наблюдений. В таблице sunspot собраны измерения, выполненные астрономами, следившими за количеством и размером пятен на поверхности Солнца.
Изучение данных¶
Загрузим необходимые библиотеки, а затем отобразим график, показывающий количество пятен на поверхности Солнца, зафиксированных между 1700 и 2000 годами.
using CSV, DataFrames
gr()
sunspot = CSV.read("sunspot.csv", DataFrame)
year = sunspot[:,1];
relNums = sunspot[:,2];
plot( year, relNums, xlabel="Год наблюдения",
ylabel="Количество пятен", title="Данные о пятнах на солнце", leg=false )
Приблизим график и попробуем более тщательно изучить циклический характер появлений солнечных пятен. Для этого построим график за первые 50 лет наблюдений.
plot( year[1:50], relNums[1:50], c=1, marker=(:circle,3),
leg=false, xlabel="Год надлюдения", ylabel="Количество пятен" )
Как мы видим, количество пятен заметно меняется, но можно предположить, что в данных есть периодический компонент, ведь пики на графике расположены с интервалом примерно в 10-15 лет.
Преобразование Фурье¶
Преобразование Фурье часто является базовым инструментом в обработке сигналов. Оно позволяет идентифицировать гармонические составляющие в данных (колебательные процессы). Используя функцию fft из предустановленной библиотеки FFTW мы преобразуем данные о пятнах на солнце из временно́го домена в частотный.
Из полученного вектора значений нужно удалить первый элемент, который хранит сумму данных. Постройте график остальной части полученного вектора. В нем содержится действительная (Re) и мнимая (Im) части комплексных коэффициентов преобразования Фурье.
using FFTW
y = fft( relNums )
y[1] = NaN
scatter( y, marker=(:o, 4, :white), markerstrokecolor=:red,
xlabel="Re(y)", ylabel="Im(y)", title="Коэффициенты Фурье", leg=false )
Коэффициенты Фурье трудно интерпретировать напрямую. Более информативным показателем вклада того или иного коэффициента в процесс является квадрат их амплитуды, который можно интерпретировать как мощность сигнала. Поскольку в векторе значений мощности каждого компонента каждое значение встречается дважды (для положительной и отрицательной шкалы частот), нам необходимо вычислить мощность любой из двух половин вектора коэффициентов. Построим спектр мощности компонентов разной частоты, показывающих количество наблюдаемых циклов за 1 год.
N = length(y);
power = abs.( y[1:Int32(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="Частота циклов (количество в году)", ylabel="Мощность", leg=false )
Вспомним что пики на графике количества пятен наблюдаются гораздо реже, чем раз в год. Чтобы график было легче читать, построим его в зависимости от периода гармонического компонента (длина цикла в годах).
График показывает, что самый заметный колебательный процесс, выражающийся в изменении количества солнечных пятен, повторяется примерно раз в 11 лет.
period = 1 ./ freq;
plot( period, power, xlimits=(0,50), xlabel="Длина цикла (годы)", ylabel="Мощность", leg=false )
Это особенно заметно после построения приближенного графика, на котором мы смогли ближе увидеть наиболее мощные составляющие сигнала.
Заключение¶
Мы изучили самый частый и простой сценарий применения быстрого преобразования Фурье (БПФ) через функцию fft библиотеки FFTW для поиска гармонических составляющих сигнала, то есть поиска в сигнале повторяющихся, циклических процессов и определения их мощности в общем сигнале.