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

Question

我有一系列的太阳黑子数据，每月计算平均太阳黑子数，我试图使用傅里叶变换将其从时间域转换到频率域。数据来自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轴的刻度。但是我知道在(11年)^-1处应该有一个峰值。

从这张图中我想了解的第二件事是为什么似乎有两条线 - 其中一条是略高于y = 0的水平线。当我将x轴边界改为：plt.xlim(0,1)时，情况更加清晰。

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

- jay

---

FFT是棘手的东西。请记住，FFT将信号分解为正弦波。仅因为您的数据每132个月有一个峰值，并不意味着您会在那里看到FFT峰值。我获取了数据，绘图显示25个月的峰值频率。顺便说一句，考虑使用np.fft.rfft，因为您有实值数据。 - Tim Roberts

谢谢您关于使用np.fft.rfft的建议！只是想问一下 - FFT可能不会每11年有一个峰值，但是numpy.fft.fftfreq应该在11年处返回一个峰值，对吗？ - jay

不，fftfreq根本不使用数据。它只是计算您的x轴值。它基本上是np.arange()除以一个常数。 - Tim Roberts

谢谢！我觉得现在我对情节和代码有了更好的理解，因为我已经稍微摆弄了一下。 - jay

Answer 1

您可以使用任何单位。请随意将采样频率表示为fs=12（每年采样次数），那么x轴的单位将为每年1个单位。或者使用fs=1（每月采样次数），单位将为每月1个单位。

您发现的额外线条来自于数据绘制方式。查看np.fft.fftfreq调用的输出。该数组的前一半包含从0到1.2e6左右的正值，另一半包含从-1.2e6到接近0的负值。通过绘制所有数据，您会得到一个从0到右边的数据线，然后是从最右边的点到最左边的点的直线，然后是其余的数据线返回到零。您的xlim调用使您无法看到绘制了一半的数据。

通常，您只需要绘制数据的前一半，只需裁剪freqs和power_spectrum数组即可。

- Cris Luengo

---

非常感谢！我现在明白这两行代码了。只有一个关于采样频率的问题，如果您不介意我问得愚蠢，那么它是fs = 12个样本/秒一年吗？还是fs = 12个样本/1年？ - jay

@jay 每月有一个样本，或者每年有12个样本。这就是你的采样频率。因此，在绘图中，你的单位将分别为1个月或1年，而不是Hz（1/s）。 - Cris Luengo

我明白这个道理，但现在我有点困惑，x轴的单位是月^-1还是年^-1，为什么总会在x=0处出现如此尖锐的峰值？因为数据中既没有每月1次也没有每年1次的基本周期性。 - jay

@jay 我已经编辑了答案，更明确地说明了采样频率和单位。现在清楚了吗？ - Cris Luengo

是的，谢谢！我的图表按预期工作了。我也相信我明白为什么在x=0处有一个强烈的峰值——有人建议我这是因为在x=0处输出是平均值。 - jay

@jay 实际上，0 频率是信号的常数部分。你可以将 0 频率设置为 0，或者在应用变换之前从信号中减去均值，结果是相同的。 - Cris Luengo

Answer 2

要将频率的步长设置为一年，确定代表此持续时间的样本数量（12），然后取参数d的倒数：

freqs = sf.rfftfreq(n=N, d=1/12)

此外，您可以通过以下几种方式改进频率比例尺：

• 使用周期而不是频率
• 反转比例尺方向，将短周期放在左边。
• 限制所绘制的频率范围（例如前50个分量），以便查看具有更长周期组件的细节。两个循环分别为11年和80-90年。

import numpy as np
import pandas as pd
import scipy.fft as sf
import matplotlib.pyplot as plt
import matplotlib.ticker as tck

# Read values, select columns, convert to numpy array
df = pd.read_excel(fp)
df = df.take([3, 1, 4], axis=1)
data = df.to_numpy()

# Sort by date, extract columns in invidual views, remove DC offset
data = data[data[:,0].argsort()]
year = data[:,1]
spots = data[:,2]
spots = spots - spots.mean()

# Get positive DFT of AQI
N = year.shape[0]
X = sf.rfft(spots) / N
freqs = sf.rfftfreq(n=N, d=1/12) # unit = 1/12 of sampling period

# Plot signal
fig, axes = plt.subplots(figsize=(15,3), ncols=2)
ax=axes[0]
ax.plot(year, spots)
ax.xaxis.set_major_locator(tck.MultipleLocator(50))
#ax.grid()

# Plot DFT
ax=axes[1]
extent = 50#N
ax.set_xlabel('period, years')
ax.stem(freqs[:extent], abs(X[:extent]))
ticks = ax.get_xticks()
ax.set_xticklabels([f'{1/tick:.1f}' if tick!=0 else '$\infty$' for tick in ticks])
ax.invert_xaxis()
ax.grid()

- mins