Engee 文档
Notebook

我们使用FFT（FFT）研究循环数据

在本演示中，我们将分析记录太阳表面斑点数量的数据，并使用快速傅立叶变换算法识别循环分量。

导言

傅里叶变换使我们能够分析数据变化的性质，当然，它们主要用于分析物理和自然过程。

我们将考虑的数据是近300年的观测历史的集合。 表中 sunspot 收集了监测太阳表面斑点数量和大小的天文学家所做的测量。

数据探索

让我们下载必要的库，然后我们将显示一个图表，显示在1700和2000之间记录的太阳表面上的斑点数量。

In [ ]: 
Pkg.add(["CSV"])
In [ ]: 
using CSV, DataFrames
gr()

sunspot = CSV.read("sunspot.csv", DataFrame)
year = sunspot[:,1];
relNums = sunspot[:,2];
plot( year, relNums, xlabel="Год наблюдения",
    ylabel="Количество пятен", title="Данные о пятнах на солнце", leg=false )
Out[0]:

让我们近似图，并尝试更彻底地研究太阳黑子出现的周期性。 为此，我们将绘制前50年观测值的图表。

In [ ]: 
plot( year[1:50], relNums[1:50], c=1, marker=(:circle,3),
    leg=false, xlabel="Год надлюдения", ylabel="Количество пятен" )
Out[0]:

正如我们所看到的，斑点的数量明显变化，但我们可以假设数据中存在周期性分量，因为图上的峰值位于大约10-15年的间隔内。

傅立叶变换

傅里叶变换通常是信号处理中的基本工具。 它允许您识别数据中的谐波分量（振荡过程）。 使用函数 fft 从预安装的库 FFTW 我们正在将太阳黑子数据从时间域转换为频率域。

从生成的值向量中，您需要删除存储数据总和的第一个元素。 绘制所得矢量的其余部分。 它包含傅立叶变换的复数系数的实部（Re）和虚部（Im）。

In [ ]: 
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 )
Out[0]:

傅立叶系数很难直接解释。 系数对该过程贡献的一个更有用的指标是其幅度的平方，可以解释为信号的功率。 由于每个值在每个分量的功率值的向量中出现两次（对于正和负频率标度），我们需要计算系数向量的两个半部分中的任何一个的功率。 让我们绘制不同频率组件的功率谱，显示在1年内观察到的周期数。

In [ ]: 
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 )
Out[0]:

回想一下，斑点数量图中的峰值观察频率远低于每年一次。 为了使图表更容易阅读，我们将根据谐波分量的周期（周期长度以年为单位）绘制它。

图表显示，最明显的振荡过程，以太阳黑子数量的变化表示，大约每11年重复一次。

In [ ]: 
period = 1 ./ freq;
plot( period, power, xlimits=(0,50), xlabel="Длина цикла (годы)", ylabel="Мощность", leg=false )
Out[0]:

在构建近似图后，这一点尤其明显，在该图中，我们能够更密切地看到信号中最强大的分量。

结论

我们通过函数研究了快速傅立叶变换（FFT）最常见和最简单的应用场景 fft 图书馆 FFTW 来搜索信号的谐波分量，即在信号中搜索重复的、循环的过程并确定它们在整体信号中的功率。