Python中的傅里叶变换时间序列答案

【问题标题】：Fourier Transform Time Series in PythonPython中的傅里叶变换时间序列
【发布时间】：2022-01-04 22:51:15
【问题描述】：

我有一个太阳黑子数量的时间序列，其中每月计算太阳黑子的平均数量，我正在尝试使用傅里叶变换从时域转换到频域。使用的数据来自https://wwwbis.sidc.be/silso/infosnmtot。我感到困惑的第一件事是如何将采样频率表示为每月一次。我是否需要将其转换为秒，例如。 1/（30 天内的秒数）？到目前为止，这是我所得到的：

fs = 1/2592000
#the sampling frequency is 1/(seconds in a month)

fourier = np.fft.fft(sn_value)
#sn_value is the mean number of sunspots measured each month
freqs = np.fft.fftfreq(sn_value.size,d=fs)

power_spectrum = np.abs(fourier)

plt.plot(freqs,power_spectrum)

plt.xlim(0,max(freqs))
plt.title("Power Spectral Density of the Sunspot Number Time Series")
plt.grid(True)

我不认为这是正确的 - 即因为我不知道 x 轴的比例是多少。但是我知道应该在 (11years)^-1 有一个高峰。

我想从这张图中的第二件事是为什么似乎有两条线 - 一条是 y=0 上方的水平线。当我将 x 轴边界更改为：plt.xlim(0,1) 时会更清楚。

我是否错误地使用了傅立叶变换函数？

【问题讨论】：

FFT 是很棘手的事情。请记住，FFT 将信号分解为正弦波。仅仅因为您的数据每 132 个月出现一次峰值并不意味着您将在那里出现 FFT 峰值。我获取了数据，该图显示了 25 个月的峰值频率。顺便说一句，您可以考虑改用np.fft.rfft，因为您有实值数据。
感谢您对使用np.fft.rfft的建议！只是问一下 - FFT 可能不会每 11 年出现一次峰值，但 numpy.fft.fftfreq 应该会在 11 年时返回一个峰值，对吗？
没有。 fftfreq 根本不使用数据。它只是计算您的 x 轴值。它基本上是np.arange() 除以一个常数。
谢谢！我想我现在理解得更好了，因为我也对情节和代码做了一些调整。

标签： python fft

【解决方案1】：

你可以使用任何你想要的单位。随意将您的采样频率表示为fs=12（样本/年），x 轴将为 1/年单位。或者使用fs=1（样本/月），则单位为1/月。

您发现的额外线条来自您绘制数据的方式。查看np.fft.fftfreq 调用的输出。该数组的前半部分包含从 0 到 1.2e6 左右的正值，另一半包含从 -1.2e6 到几乎为 0 的负值。通过绘制所有数据，您会得到一条从 0 到右侧的数据线，然后是从最右边点到最左边点的直线，然后将其余数据线归零。您的xlim 调用成功了，因此您看不到绘制的一半数据。

通常您只会绘制数据的前半部分，只需裁剪 freqs 和 power_spectrum 数组。

【讨论】：

非常感谢！我现在明白这两条线了。如果您不介意我很愚蠢，只是一个关于采样频率的快速问题，但是一年中 fs = 12 个样本/秒吗？还是 fs = 12samples/1year？
@jay 您每个月有一个样本，或者每年有 12 个样本。那是您的采样频率。因此，图中的单位将分别为 1/月或 1/年，而不是 Hz (1/s)。
我认为这是有道理的——虽然我现在有点困惑，如果 x 轴的单位是月^-1 或年^-1，为什么总是会在x=0？因为 1/月和 1/年都不应该是数据中的基本周期。
@jay 我已经编辑了答案，以便更明确地了解采样频率和单位。现在清楚了吗？
@jay 实际上，0 频率是信号的恒定部分。您可以将频率设置为 0 到 0，也可以在应用变换之前从信号中减去均值，结果是一样的。