用Python绘制快速傅里叶变换

Question

我可以访问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的参数改为我的数据集并绘制出来时，我得到了非常奇怪的结果，看起来频率的缩放可能是关闭的。我不确定。

这是我试图快速傅里叶变换的数据的一个粘贴箱

http://pastebin.com/0WhjjMkbhttp://pastebin.com/ksM4FvZS

当我对整件事做快速傅里叶变换时，它只在0处有一个巨大的尖峰，其他的都没有

这是我的代码：

## 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]

Answer 1

您的高峰值是由于您的信号的直流（非变化，即频率=0）部分。这是一个规模问题。如果要查看非直流频率内容，为了可视化，可能需要从偏移量1而不是从信号FFT的偏移量0绘制。

修改上面@PaulH给出的示例

import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack

# Number of samplepoints
N = 600
# sample spacing
T = 1.0 / 800.0
x = np.linspace(0.0, N*T, N)
y = 10 + np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x)
yf = scipy.fftpack.fft(y)
xf = np.linspace(0.0, 1.0/(2.0*T), N/2)

plt.subplot(2, 1, 1)
plt.plot(xf, 2.0/N * np.abs(yf[0:N/2]))
plt.subplot(2, 1, 2)
plt.plot(xf[1:], 2.0/N * np.abs(yf[0:N/2])[1:])

输出图：

另一种方法是，以对数比例可视化数据：

使用：

plt.semilogy(xf, 2.0/N * np.abs(yf[0:N/2]))

将显示：

Answer 2

所以我在一个IPython笔记本上运行了一个功能相当的代码：

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack

# 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 = scipy.fftpack.fft(y)
xf = np.linspace(0.0, 1.0/(2.0*T), N/2)

fig, ax = plt.subplots()
ax.plot(xf, 2.0/N * np.abs(yf[:N//2]))
plt.show()

我得到了我认为非常合理的结果。

自从我在工科学校考虑信号处理以来，这比我想承认的要长，但50和80的峰值正是我所期望的。那有什么问题？

针对发布的原始数据和评论

这里的问题是你没有周期性的数据。您应该始终检查输入任何算法的数据，以确保它是适当的。

import pandas
import matplotlib.pyplot as plt
#import seaborn
%matplotlib inline

# the OP's data
x = pandas.read_csv('http://pastebin.com/raw.php?i=ksM4FvZS', skiprows=2, header=None).values
y = pandas.read_csv('http://pastebin.com/raw.php?i=0WhjjMkb', skiprows=2, header=None).values
fig, ax = plt.subplots()
ax.plot(x, y)

Answer 3

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函数

Y    = numpy.fft.fft(y)
freq = numpy.fft.fftfreq(len(y), t[1] - t[0])

pylab.figure()
pylab.plot( freq, numpy.abs(Y) )
pylab.figure()
pylab.plot(freq, numpy.angle(Y) )
pylab.show()

这应该能解决你的问题。