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 轴的比例是多少。但是我知道在 (11years)^-1 应该有一个高峰。

第二个我想知道的事情是为什么看起来有两条线 - 一条是 y=0 上方的水平线。当我将 x 轴边界更改为：plt.xlim(0,1).

时更清楚

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

Answer 1

您可以使用任何您想要的单位。随意将您的采样频率表示为

fs=12

（样本/年），x 轴将以 1/年为单位。或使用

fs=1

（样本/月），则单位为1/月。

您发现的额外线条来自您绘制数据的方式。查看

np.fft.fftfreq

调用的输出。该数组的前半部分包含从 0 到 1.2e6 左右的正值，另一半包含从 -1.2e6 到几乎为 0 的负值。通过绘制所有数据，您会得到一条从 0 到右边的数据线，然后是从最右边的点到最左边的点直线，然后数据线的其余部分回零。你的

xlim

电话成功了，所以你看不到一半的数据绘制。

通常你只会绘制数据的前半部分，只需裁剪

freqs

和

power_spectrum

数组。

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()