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)

---

第二次编辑：

侦探的回答让我产生了我正在寻找的情节。对任何感兴趣的人来说，情节如下。

未经过滤的数据系列-

在下面，对上述数据序列执行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

^{}函数返回完整的N点谱（对于实值输入，包括多余的谱上半部分），而您的频率轴xf被构造为只包含N//2点的谱下半部分。您的错误与这些xf和yf数组大小之间的不匹配有关。由于冗余，您可以使用yf[0:N//2]排除yf中频谱的上半部分。

还要注意数组yf包含复数。要显示光谱图输出，应采用绝对值：

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

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