Python：使用 SciPy 文档对 .csv 值执行 FFT
Python: Performing FFT on .csv values using SciPy documentation
我想对数据序列进行快速傅里叶变换。该系列包含每日地震振幅的值，在 407 天内连续采样。我想在这个数据集中搜索任何周期性周期。
我最初尝试使用 SciPy 文档：https://docs.scipy.org/doc/scipy/reference/tutorial/fftpack.html. Similar to this question (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 生成的图 -
我希望这是正确的输出，并且我已经正确标记了 axes/understood 第二个图代表的内容。我还想知道为什么傅里叶谱的上半部分是多余的。
作为参考，我的完整（更正）代码如下：
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()
fft 函数 returns 完整 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（或cycles/day）作为频率，在这种情况下您可以使用T = 1.
import pandas as pd
import numpy as np
from numpy.fft import rfft, rfftfreq
import matplotlib.pyplot as plt
t=pd.read_csv('C:\Users\trial\Desktop\EW.csv',usecols=[0])
a=pd.read_csv('C:\Users\trial\Desktop\EW.csv',usecols=[1])
n=len(a)
dt=0.02 #time increment in each data
acc=a.values.flatten() #to convert DataFrame to 1D array
#acc value must be in numpy array format for half way mirror calculation
fft=rfft(acc)*dt
freq=rfftfreq(n,d=dt)
FFT=abs(fft)
plt.plot(freq,FFT)
我想对数据序列进行快速傅里叶变换。该系列包含每日地震振幅的值，在 407 天内连续采样。我想在这个数据集中搜索任何周期性周期。
我最初尝试使用 SciPy 文档：https://docs.scipy.org/doc/scipy/reference/tutorial/fftpack.html. Similar to this question (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 生成的图 -
我希望这是正确的输出，并且我已经正确标记了 axes/understood 第二个图代表的内容。我还想知道为什么傅里叶谱的上半部分是多余的。
作为参考，我的完整（更正）代码如下：
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()
fft 函数 returns 完整 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（或cycles/day）作为频率，在这种情况下您可以使用T = 1.
import pandas as pd
import numpy as np
from numpy.fft import rfft, rfftfreq
import matplotlib.pyplot as plt
t=pd.read_csv('C:\Users\trial\Desktop\EW.csv',usecols=[0])
a=pd.read_csv('C:\Users\trial\Desktop\EW.csv',usecols=[1])
n=len(a)
dt=0.02 #time increment in each data
acc=a.values.flatten() #to convert DataFrame to 1D array
#acc value must be in numpy array format for half way mirror calculation
fft=rfft(acc)*dt
freq=rfftfreq(n,d=dt)
FFT=abs(fft)
plt.plot(freq,FFT)