使用 FFT（FFT）探索循环数据¶

在本演示中，我们使用快速傅立叶变换算法分析包含太阳表面观测到的太阳黑子数量年表的数据，并识别周期成分。

简介¶

傅立叶变换可以让我们分析数据的变化模式，当然，最早的应用是分析物理和自然过程。

我们要查看的数据是近 300 年观测历史的汇编。表sunspot 总结了天文学家对太阳表面太阳黑子的数量和大小进行的测量。

查看数据¶

让我们加载必要的库，然后显示一张 1700 年至 2000 年太阳表面记录的黑子数量图。

Pkg.add(["CSV"])

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 年。

傅立叶变换¶

傅立叶变换通常是信号处理的基本工具。它可以识别数据中的谐波成分（振荡过程）。利用预装库FFTW 中的函数fft ，我们将太阳黑子数据从时域转换到频域。

我们需要从得到的数值矢量中移除第一个元素，它存储了数据的总和。将生成向量的其余部分绘制成图。它包含傅立叶变换复数系数的实（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 )

这一点在绘制近似图后尤为明显，我们可以在近似图上更近距离地看到信号中最强烈的成分。

结论¶

我们研究了通过FFTW 库的fft 函数应用快速傅立叶变换（FFT）来搜索信号谐波成分的最常见和最简单的情况，即查找信号中的重复和循环过程，并确定它们在整个信号中的功率。