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

Question

我想对数据序列执行快速傅立叶变换。该系列包含每日地震振幅的值，在 407 天内连续采样。我想在此数据集中搜索任何周期性周期。

我在这里使用 SciPy 文档进行了初步尝试:https://docs.scipy.org/doc/scipy/reference/tutorial/fftpack.html .与这个问题 (link) 类似，然后我将 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[0:N//2] 排除 yf 中频谱的上半部分。

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

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

最后就采样周期而言，如果您要使用秒作为采样周期并使用 Hz 作为频率，则您应该使用 T = 86400(因为您每 1 天或86400 秒)。您也可以选择使用天作为采样周期，使用天^-1(或周期/天)作为频率，在这种情况下您可以使用T = 1。