在 Python 中绘制快速傅立叶变换
声明：本页面是StackOverFlow热门问题的中英对照翻译，遵循CC BY-SA 4.0协议，如果您需要使用它，必须同样遵循CC BY-SA许可，注明原文地址和作者信息，同时你必须将它归于原作者（不是我）：StackOverFlow
原文地址: http://stackoverflow.com/questions/25735153/
Warning: these are provided under cc-by-sa 4.0 license. You are free to use/share it, But you must attribute it to the original authors (not me):
StackOverFlow
Plotting a Fast Fourier Transform in Python
提问by user3123955
I have access to numpy and scipy and want to create a simple FFT of a dataset. I have two lists one that is y values and the other is timestamps for those y values.
我可以访问 numpy 和 scipy 并想创建一个简单的数据集 FFT。我有两个列表，一个是 y 值，另一个是那些 y 值的时间戳。
What is the simplest way to feed these lists into a scipy or numpy method and plot the resulting FFT?
将这些列表输入 scipy 或 numpy 方法并绘制生成的 FFT 的最简单方法是什么？
I have looked up examples, but they all rely on creating a set of fake data with some certain number of data points, and frequency, etc. and doesn't really show how to do it with just a set of data and the corresponding timestamps.
我查了一些例子，但它们都依赖于创建一组具有一定数量的数据点和频率等的假数据，并没有真正展示如何只用一组数据和相应的时间戳来做到这一点.
I have tried the following example:
我尝试了以下示例：
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()
But when i change the argument of fft to my data set and plot it i get extremely odd results, it appears the scaling for the frequency may be off. i am unsure.
但是当我将 fft 的参数更改为我的数据集并绘制它时，我得到了非常奇怪的结果，似乎频率的缩放可能关闭。我不确定。
Here is a pastebin of the data i am attempting to FFT
这是我尝试 FFT 的数据的粘贴箱
http://pastebin.com/0WhjjMkbhttp://pastebin.com/ksM4FvZS
http://pastebin.com/0WhjjMkb http://pastebin.com/ksM4FvZS
When i do an fft on the whole thing it just has a huge spike at zero and nothing else
当我对整件事做一个 fft 时，它只是在零处有一个巨大的尖峰，没有别的
Here is my code:
这是我的代码：
## 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)')
spacing is just equal to xInterp[1]-xInterp[0]
间距刚好等于 xInterp[1]-xInterp[0]
采纳答案by Paul H
So I run a functionally equivalent form of your code in an IPython notebook:
因此，我在 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()
I get what I believe to be very reasonable output.
我得到了我认为非常合理的输出。
It's been longer than I care to admit since I was in engineering school thinking about signal processing, but spikes at 50 and 80 are exactly what I would expect. So what's the issue?
自从我在工程学校考虑信号处理以来，我承认的时间比我想承认的要长，但是 50 和 80 的峰值正是我所期望的。那么问题是什么？
In response to the raw data and comments being posted
回应发布的原始数据和评论
The problem here is that you don't have periodic data. You should always inspect the data that you feed into anyalgorithm to make sure that it's appropriate.
这里的问题是您没有定期数据。您应该始终检查您输入任何算法的数据，以确保它是合适的。
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)
回答by ssm
The important thing about fft is that it can only be applied to data in which the timestamp is uniform (i.e.uniform sampling in time, like what you have shown above).
关于 fft 的重要一点是它只能应用于时间戳是统一的数据（即时间上的统一采样，就像你上面展示的那样）。
In case of non-uniform sampling, please use a function for fitting the data. There are several tutorials and functions to choose from:
如果采样不均匀，请使用拟合数据的函数。有几个教程和功能可供选择：
https://github.com/tiagopereira/python_tips/wiki/Scipy%3A-curve-fittinghttp://docs.scipy.org/doc/numpy/reference/generated/numpy.polyfit.html
https://github.com/tiagopereira/python_tips/wiki/Scipy%3A-curve-fit http://docs.scipy.org/doc/numpy/reference/generated/numpy.polyfit.html
If fitting is not an option, you can directly use some form of interpolation to interpolate data to a uniform sampling:
如果拟合不是一个选项，您可以直接使用某种形式的插值将数据插值到统一采样：
https://docs.scipy.org/doc/scipy-0.14.0/reference/tutorial/interpolate.html
https://docs.scipy.org/doc/scipy-0.14.0/reference/tutorial/interpolate.html
When you have uniform samples, you will only have to wory about the time delta (t[1] - t[0]) of your samples. In this case, you can directly use the fft functions
当您有统一的样本时，您只需担心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()
This should solve your problem.
这应该可以解决您的问题。
回答by hesham_EE
The high spike that you have is due to the DC (non-varying, i.e. freq = 0) portion of your signal. It's an issue of scale. If you want to see non-DC frequency content, for visualization, you may need to plot from the offset 1 not from offset 0 of the FFT of the signal.
您拥有的高尖峰是由于信号的 DC（不变，即 freq = 0）部分造成的。这是规模问题。如果您想查看非 DC 频率内容，为了可视化，您可能需要从信号 FFT 的偏移 1 而非偏移 0 进行绘图。
Modifying the example given above by @PaulH
修改@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:])
The output plots:
输出图：
Another way, is to visualize the data in log scale:
另一种方法是以对数比例可视化数据：
Using:
使用：
plt.semilogy(xf, 2.0/N * np.abs(yf[0:N/2]))
Will show:
将会呈现：
回答by ewerlopes
Just as a complement to the answers already given, I would like to point out that often it is important to play with the size of the bins for the FFT. It would make sense to test a bunch of values and pick the one that makes more sense to your application. Often, it is in the same magnitude of the number of samples. This was as assumed by most of the answers given, and produces great and reasonable results. In case one wants to explore that, here is my code version:
作为对已经给出的答案的补充，我想指出，通常使用 FFT 的 bin 大小很重要。测试一堆值并选择对您的应用程序更有意义的值是有意义的。通常，它与样本数量的数量级相同。这是大多数给出的答案所假设的，并产生了很好且合理的结果。如果有人想探索这一点，这是我的代码版本：
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack
fig = plt.figure(figsize=[14,4])
N = 600 # Number of samplepoints
Fs = 800.0
T = 1.0 / Fs # N_samps*T (#samples x sample period) is the sample spacing.
N_fft = 80 # Number of bins (chooses granularity)
x = np.linspace(0, N*T, N) # the interval
y = np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x) # the signal
# removing the mean of the signal
mean_removed = np.ones_like(y)*np.mean(y)
y = y - mean_removed
# Compute the fft.
yf = scipy.fftpack.fft(y,n=N_fft)
xf = np.arange(0,Fs,Fs/N_fft)
##### Plot the fft #####
ax = plt.subplot(121)
pt, = ax.plot(xf,np.abs(yf), lw=2.0, c='b')
p = plt.Rectangle((Fs/2, 0), Fs/2, ax.get_ylim()[1], facecolor="grey", fill=True, alpha=0.75, hatch="/", zorder=3)
ax.add_patch(p)
ax.set_xlim((ax.get_xlim()[0],Fs))
ax.set_title('FFT', fontsize= 16, fontweight="bold")
ax.set_ylabel('FFT magnitude (power)')
ax.set_xlabel('Frequency (Hz)')
plt.legend((p,), ('mirrowed',))
ax.grid()
##### Close up on the graph of fft#######
# This is the same histogram above, but truncated at the max frequence + an offset.
offset = 1 # just to help the visualization. Nothing important.
ax2 = fig.add_subplot(122)
ax2.plot(xf,np.abs(yf), lw=2.0, c='b')
ax2.set_xticks(xf)
ax2.set_xlim(-1,int(Fs/6)+offset)
ax2.set_title('FFT close-up', fontsize= 16, fontweight="bold")
ax2.set_ylabel('FFT magnitude (power) - log')
ax2.set_xlabel('Frequency (Hz)')
ax2.hold(True)
ax2.grid()
plt.yscale('log')
the output plots:
输出图：
回答by Foad
There are already great solutions on this page, but all have assumed the dataset is uniformly/evenly sampled/distributed. I will try to provide a more general example of randomly sampled data. I will also use this MATLAB tutorialas an example:
这个页面上已经有很好的解决方案，但都假设数据集是均匀/均匀采样/分布的。我将尝试提供一个更一般的随机采样数据示例。我还将使用此 MATLAB 教程作为示例：
Adding the required modules:
添加所需的模块：
import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack
import scipy.signal
Generating sample data:
生成样本数据：
N = 600 # number of samples
t = np.random.uniform(0.0, 1.0, N) # assuming the time start is 0.0 and time end is 1.0
S = 1.0 * np.sin(50.0 * 2 * np.pi * t) + 0.5 * np.sin(80.0 * 2 * np.pi * t)
X = S + 0.01 * np.random.randn(N) # adding noise
Sorting the data set:
对数据集进行排序：
order = np.argsort(t)
ts = np.array(t)[order]
Xs = np.array(X)[order]
Resampling:
重采样：
T = (t.max() - t.min()) / N # average period
Fs = 1 / T # average sample rate frequency
f = Fs * np.arange(0, N // 2 + 1) / N; # resampled frequency vector
X_new, t_new = scipy.signal.resample(Xs, N, ts)
plotting the data and resampled data:
绘制数据和重新采样的数据：
plt.xlim(0, 0.1)
plt.plot(t_new, X_new, label="resampled")
plt.plot(ts, Xs, label="org")
plt.legend()
plt.ylabel("X")
plt.xlabel("t")
now calculating the fft:
现在计算fft：
Y = scipy.fftpack.fft(X_new)
P2 = np.abs(Y / N)
P1 = P2[0 : N // 2 + 1]
P1[1 : -2] = 2 * P1[1 : -2]
plt.ylabel("Y")
plt.xlabel("f")
plt.plot(f, P1)
P.S.I finally got time to implement a more canonical algorithm to get a Fourier transform of unevenly distributed data. You may see the code, description, and example Jupyter notebook here.
PS我终于有时间实现一个更规范的算法来获得不均匀分布数据的傅立叶变换。您可以在此处查看代码、描述和示例 Jupyter 笔记本。
回答by YoniChechik
I've built a function that deals with plotting FFT of real signals. The extra bonus in my function relative to the messages above is that you get the ACTUAL amplitude of the signal. Also, because of the assumption of a real signal, the FFT is symmetric so we can plot only the positive side of the x axis:
我已经构建了一个函数来处理绘制真实信号的 FFT。相对于上述消息，我的函数中的额外好处是您可以获得信号的实际幅度。此外，由于假设是真实信号，FFT 是对称的，因此我们只能绘制 x 轴的正侧：
import matplotlib.pyplot as plt
import numpy as np
import warnings
def fftPlot(sig, dt=None, plot=True):
# here it's assumes analytic signal (real signal...)- so only half of the axis is required
if dt is None:
dt = 1
t = np.arange(0, sig.shape[-1])
xLabel = 'samples'
else:
t = np.arange(0, sig.shape[-1]) * dt
xLabel = 'freq [Hz]'
if sig.shape[0] % 2 != 0:
warnings.warn("signal prefered to be even in size, autoFixing it...")
t = t[0:-1]
sig = sig[0:-1]
sigFFT = np.fft.fft(sig) / t.shape[0] # divided by size t for coherent magnitude
freq = np.fft.fftfreq(t.shape[0], d=dt)
# plot analytic signal - right half of freq axis needed only...
firstNegInd = np.argmax(freq < 0)
freqAxisPos = freq[0:firstNegInd]
sigFFTPos = 2 * sigFFT[0:firstNegInd] # *2 because of magnitude of analytic signal
if plot:
plt.figure()
plt.plot(freqAxisPos, np.abs(sigFFTPos))
plt.xlabel(xLabel)
plt.ylabel('mag')
plt.title('Analytic FFT plot')
plt.show()
return sigFFTPos, freqAxisPos
if __name__ == "__main__":
dt = 1 / 1000
# build a signal within nyquist - the result will be the positive FFT with actual magnitude
f0 = 200 # [Hz]
t = np.arange(0, 1 + dt, dt)
sig = 1 * np.sin(2 * np.pi * f0 * t) + \
10 * np.sin(2 * np.pi * f0 / 2 * t) + \
3 * np.sin(2 * np.pi * f0 / 4 * t) +\
7.5 * np.sin(2 * np.pi * f0 / 5 * t)
# res in freqs
fftPlot(sig, dt=dt)
# res in samples (if freqs axis is unknown)
fftPlot(sig)
