Python：使用 SciPy 文档对 .csv 值执行 FFT

Question

我想对数据系列执行快速傅里叶变换。该系列包含 407 天内一致采样的每日地震振幅值。我想在这个数据集中搜索任何周期性周期。

我在这里使用 SciPy 文档进行了初步尝试：https://docs.scipy.org/doc/scipy/reference/tutorial/fftpack.html。类似于这个问题（链接），然后我将 y 的参数从人工正弦函数更改为我的数据集。

但是，我收到以下错误：

ValueError: x and y must have same first dimension, but have shapes (203,) and (407, 1)

我希望能帮助您理解为什么会出现此错误以及如何修复它。

我还希望获得有关 FFT 处理数据集所需的正确频率和采样输入值的帮助。我的数据集中有 407 个值，每个值代表一天。因此，我定义了 N（样本点数）= 407，T（样本间距）= 1 / 84600（1 / 一天的秒数）。那是对的吗？

这是我的完整代码：

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft, ifft
import pandas as pd

# Import csv file
df = pd.read_csv('rsam_2016-17_fft_test.csv', index_col=['DateTime'], parse_dates=['DateTime'])
print(df.head())

#plot data
plt.figure(figsize=(12,4))
df.plot(linestyle = '', marker = '*', color='r')
plt.show()

#FFT
#number of sample points
N = 407
#frequency of signal
T = 1 / 84600
#create x-axis for time length of signal
x = np.linspace(0, N*T, N)
#create array that corresponds to values in signal
y = df
#perform FFT on signal
yf = fft(y)
#create new x-axis: frequency from signal
xf = np.linspace(0.0, 1.0/(2.0*T), N//2)
#plot results
plt.plot(xf, yf)
plt.grid()
plt.show()

非常感谢您的帮助！

---

编辑：完整的错误如下：

Traceback (most recent call last):

File "<ipython-input-40-c090e0039ba9>", line 1, in <module>
runfile('/Users/an16975/Desktop/PhD/Code/FFT/fft_test.py', wdir='/Users/an16975/Desktop/PhD/Code/FFT')

File "/Users/an16975/anaconda3/lib/python3.6/site-packages/spyder/utils/site/sitecustomize.py", line 710, in runfile
execfile(filename, namespace)

File "/Users/an16975/anaconda3/lib/python3.6/site-packages/spyder/utils/site/sitecustomize.py", line 101, in execfile
exec(compile(f.read(), filename, 'exec'), namespace)

File "/Users/an16975/Desktop/PhD/Code/FFT/fft_test.py", line 36, in <module>
plt.plot(xf, yf)

File "/Users/an16975/anaconda3/lib/python3.6/site-packages/matplotlib/pyplot.py", line 3317, in plot
ret = ax.plot(*args, **kwargs)

File "/Users/an16975/anaconda3/lib/python3.6/site-packages/matplotlib/__init__.py", line 1898, in inner
return func(ax, *args, **kwargs)

File "/Users/an16975/anaconda3/lib/python3.6/site-packages/matplotlib/axes/_axes.py", line 1406, in plot
for line in self._get_lines(*args, **kwargs):

File "/Users/an16975/anaconda3/lib/python3.6/site-packages/matplotlib/axes/_base.py", line 407, in _grab_next_args
for seg in self._plot_args(remaining, kwargs):

File "/Users/an16975/anaconda3/lib/python3.6/site-packages/matplotlib/axes/_base.py", line 385, in _plot_args
x, y = self._xy_from_xy(x, y)

File "/Users/an16975/anaconda3/lib/python3.6/site-packages/matplotlib/axes/_base.py", line 244, in _xy_from_xy
"have shapes {} and {}".format(x.shape, y.shape))

ValueError: x and y must have same first dimension, but have shapes (203,) and (407, 1)

---

第二次编辑：

SleuthEye 的答案使我能够生成我正在寻找的图。对于任何有兴趣的人，情节如下。

未过滤的数据系列 -

下面是对上述数据系列执行 FFT 生成的图 -

我希望这是正确的输出，并且我已经正确标记了轴/理解第二个图代表的内容。我还想知道为什么傅里叶谱的上部是多余的。

作为参考，我的完整（更正后）代码如下：

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft, ifft
import pandas as pd

# Import csv file
df = pd.read_csv('rsam_2016-17_fft_test.csv', index_col=['DateTime'], parse_dates=['DateTime'])
print(df.head())

#plot data
plt.figure(figsize=(12,4))
df.plot(linestyle = '', marker = '*', color='r')
plt.savefig('rsam_2016_2017_snippetforfft.jpg')
plt.show()

#FFT
#number of sample points
N = 407
#frequency of signal (in days)
T = 1
#create x-axis for time length of signal
x = np.linspace(0, N*T, N)
#create array that corresponds to values in signal
y = df
#perform FFT on signal
yf = fft(y)
#create new x-axis: frequency from signal
xf = np.linspace(0.0, 1.0/(2.0*T), N//2)
#plot results
plt.plot(xf, yf[0:N//2], label = 'signal')
plt.grid()
plt.xlabel('Frequency (days)')
plt.ylabel(r'Spectral Amplitude')
plt.legend(loc=1)
plt.savefig('rsam_2016_2017_snippet_fft_firstresult.jpg')
plt.show()

Answer 1

该fft函数返回完整的N点频谱（对于实值输入，包括频谱的冗余上半部分），而频率轴的xf构造是仅覆盖频谱的下半部N//2分。xf您的错误与这些大小和数组大小之间的不匹配有关yf。由于冗余，您可以在yf使用时排除频谱的上半部分yf[0:N//2]。

另请注意，该数组yf包含复数。要显示频谱图输出，您应该取绝对值：

plt.plot(xf, abs(yf[0:N//2]))

最后，就采样周期而言，如果您要使用秒作为采样周期，使用赫兹作为频率，则应该使用T = 86400（因为每 1 天或 86400 秒就有一个数据点）。您还可以使用天作为采样周期，使用天^-1（或周期/天）作为频率，在这种情况下，您将使用T = 1。