我可以访问NumPy和SciPy,并希望创建一个数据集的简单FFT。我有两个列表,一个是y
值,另一个是这些y
值的时间戳。
将这些列表输入SciPy或NumPy方法并绘制FFT结果的最简单方法是什么?
我查过一些例子,但它们都依赖于创建一组具有一定数量的数据点和频率等的假数据,而没有真正展示如何仅使用一组数据和相应的时间戳来做到这一点。
我试过下面的例子:
from scipy.fftpack import fft
# Number of samplepoints
N = 600
# Sample spacing
T = 1.0 / 800.0
x = np.linspace(0.0, N*T, N)
y = np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x)
yf = fft(y)
xf = np.linspace(0.0, 1.0/(2.0*T), N/2)
import matplotlib.pyplot as plt
plt.plot(xf, 2.0/N * np.abs(yf[0:N/2]))
plt.grid()
plt.show()
但是,当我将fft
的参数更改为我的数据集并绘制它时,我得到了非常奇怪的结果,而且似乎频率的缩放可能是错误的。
这是我尝试进行FFT的数据的粘贴框
http://pastebin.com/0WhjjMkbhttp://pastebin.com/ksM4FvZS显示器
当我在整个过程中使用fft()
时,它只是在零处有一个巨大的尖峰,其他什么都没有。
下面是我的代码:
## Perform FFT with SciPy
signalFFT = fft(yInterp)
## Get power spectral density
signalPSD = np.abs(signalFFT)**2
## Get frequencies corresponding to signal PSD
fftFreq = fftfreq(len(signalPSD), spacing)
## Get positive half of frequencies
i = fftfreq>0
##
plt.figurefigsize = (8, 4)
plt.plot(fftFreq[i], 10*np.log10(signalPSD[i]));
# plt.xlim(0, 100);
plt.xlabel('Frequency [Hz]');
plt.ylabel('PSD [dB]')
间距正好等于xInterp[1]-xInterp[0]
。
7条答案
按热度按时间mrphzbgm1#
因此,我在IPython笔记本中运行了一个与您的代码在功能上等效的形式:
我得到了我认为非常合理的输出。
自从我上工程学院开始思考信号处理问题以来,我已经很久没有这样做了,但50和80的峰值正是我所期望的。
作为对发布的原始数据和评论的回应
这里的问题是你没有周期性的数据,你应该总是检查你输入到 * 任何 * 算法中的数据,以确保它是合适的。
siotufzp2#
关于fft的重要一点是,它只能应用于时间戳一致的数据(* 即 * 时间上的一致采样,如上面所示)。
如果是非均匀采样,请使用函数来拟合数据。有几个教程和函数可供选择:
https://github.com/tiagopereira/python_tips/wiki/Scipy%3A-curve-fittinghttp://docs.scipy.org/doc/numpy/reference/generated/numpy.polyfit.html显示器
如果拟合不是一个选项,则可以直接使用某种形式的插值将数据插值为均匀采样:
https://docs.scipy.org/doc/scipy-0.14.0/reference/tutorial/interpolate.html
当你有均匀的样本时,你只需要担心你的样本的时间增量(
t[1] - t[0]
),在这种情况下,你可以直接使用fft函数这应该能解决你的问题。
t30tvxxf3#
高尖峰是由于信号的DC(不变,即freq = 0)部分造成的。这是一个标度问题。如果要查看非DC频率内容,为了可视化,您可能需要从信号FFT的偏移1而不是偏移0绘图。
修改上面@PaulH给出的示例
输出出图:
另一种方法是以对数刻度可视化数据:
使用:
将显示:
am46iovg4#
我已经建立了一个函数来处理真实的信号的FFT绘图。相对于前面的答案,我的函数的额外好处是你可以得到信号的 * 实际 * 幅度。
此外,由于假设为真实的信号,FFT是对称的,因此我们只能绘制x轴的正侧:
sauutmhj5#
作为对上述答案的补充,我想指出的是,通常情况下,调整FFT的bin大小很重要。测试一组值,然后选择对应用更有意义的值,这是很有意义的。通常情况下,它与样本数的大小相同。大多数答案都是这样假设的。并产生了很好的合理的结果。如果有人想探索这一点,这里是我的代码版本:
输出曲线:
xqnpmsa86#
我写这个额外的答案是为了解释使用FFT时尖峰扩散的起源,特别是讨论我在某些地方不同意的scipy.fftpack教程。
在本例中,记录时间为
tmax=N*T=0.75
。信号为sin(50*2*pi*x) + 0.5*sin(80*2*pi*x)
。频率信号应包含频率为50
和80
的两个尖峰,振幅为1
和0.5
。但是,如果分析的信号不具有整数个周期,则由于信号截断,可能出现扩散:50*tmax=37.5
=〉频率50
不是1/tmax
=〉的倍数,由于在该频率处的信号截断,存在扩散。80*tmax=60
=〉频率80
是1/tmax
=〉的倍数,由于在该频率处的信号截断,因此没有扩散。以下代码分析的信号与教程(
sin(50*2*pi*x) + 0.5*sin(80*2*pi*x)
)中的信号相同,但略有不同:1.原始的scipy.fftpack示例。
1.原始scipy.fftpack示例具有整数个信号周期(
tmax=1.0
代替0.75
以避免截断扩散)。1.原始scipy.fftpack示例具有整数个信号周期,其中日期和频率取自FFT理论。
代码:
输出量:
在这里,即使使用整数个周期,仍然会存在一些扩散。这种行为是由于scipy.fftpack教程中日期和频率的错误定位所致。因此,在离散傅里叶变换理论中:
t=0,T,...,(N-1)*T
进行评估,其中T为采样周期,信号的总持续时间为tmax=N*T
。注意,我们在tmax-T
处停止。f=0,df,...,(N-1)*df
,其中df=1/tmax=1/(N*T)
为采样频率。信号的所有谐波均应为采样频率的倍数,以避免扩散。在上面的例子中,你可以看到使用
arange
代替linspace
可以避免频谱中的额外扩散。使用linspace
形式还导致尖峰的偏移,所述尖峰位于比它们应该在的频率稍高的频率处,如在第一幅图中可以看到的,其中尖峰在频率为50
和80
。我将得出结论,应该用下面的代码替换用法示例(在我看来,这不容易引起误解):
输出(第二个尖峰不再扩散):
我认为这个答案仍然带来了一些关于如何正确应用离散傅立叶变换的额外解释。显然,我的答案太长了,而且总是有额外的东西要说(例如,ewerlopes简要地谈到了aliasing,还有很多关于窗口的内容),所以我就不多说了。
我认为在应用离散傅里叶变换时,深入理解它的原理是非常重要的,因为我们都知道,很多人在应用它时,为了得到他们想要的东西,到处添加因子。
mrwjdhj37#
本页中已经有了很好的解决方案,但所有的解决方案都假设数据集是均匀/均匀采样/分布的。我将尝试提供一个更一般的随机采样数据的示例。我还将使用this MATLAB tutorial作为示例:
添加所需模块:
正在生成示例数据:
排序数据集:
重新取样:
绘制数据和重采样数据:
现在计算FFT:
**P.S.**我终于有时间实现一个更规范的算法来获得不均匀分布数据的傅立叶变换。您可以查看代码、描述和示例Jupyter笔记本here。