在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的参数更改为我的数据集并绘制它时,我得到了非常奇怪的结果,看起来频率的缩放可能会关闭.我不确定.

这是我试图FFT的数据的pastebin

http://pastebin.com/0WhjjMkb http://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]

Answer 1

所以我在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 2

关于fft的重要之处在于它只能应用于时间戳统一的数据(即时间上的均匀采样,就像上面显示的那样).

如果采样不均匀,请使用函数来拟合数据.有几个教程和功能可供选择:

https://github.com/tiagopereira/python_tips/wiki/Scipy%3A-curve-fitting http://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]样本的时间delta().在这种情况下,您可以直接使用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()

这应该可以解决您的问题.

Answer 3

您的高峰值是由于信号的DC(非变化,即freq = 0)部分.这是一个规模问题.如果要查看非DC频率内容,可视化,您可能需要从偏移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 4

我已经构建了一个函数来处理绘制真实信号的 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 preferred 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 frequence 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)

    # Result in frequencies
    fftPlot(sig, dt=dt)
    # Result in samples (if the frequencies axis is unknown)
    fftPlot(sig)

Answer 5

作为对已经给出的答案的补充，我想指出的是，经常需要考虑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')

输出图：

Answer 6

这个页面上已经有很好的解决方案，但都假设数据集是均匀/均匀采样/分布的。我将尝试提供一个更一般的随机采样数据示例。我还将使用此 MATLAB 教程作为示例：

添加所需的模块：

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

生成样本数据：

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

对数据集进行排序：

order = np.argsort(t)
ts = np.array(t)[order]
Xs = np.array(X)[order]

重采样：

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)

绘制数据和重采样数据：

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")

现在计算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)

PS我终于有时间实现一个更规范的算法来获得不均匀分布数据的傅立叶变换。您可以在此处查看代码、描述和示例 Jupyter 笔记本。

Answer 7

我写了这个额外的答案来解释使用 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。但是，如果分析的信号没有整数个周期，则由于信号的截断可能会出现扩散：

• Pike 1：50*tmax=37.5=> 频率50不是1/tmax=>由于该频率的信号截断而导致的扩散的倍数。
• 派克 2：80*tmax=60=> 频率80是1/tmax=>由于该频率的信号截断而没有扩散的倍数。

这是一个分析与教程 ( sin(50*2*pi*x) + 0.5*sin(80*2*pi*x)) 中相同的信号的代码，但略有不同：

1. 原始 scipy.fftpack 示例。
2. 带有整数个信号周期的原始 scipy.fftpack 示例（tmax=1.0而不是0.75为了避免截断扩散）。
3. 原始 scipy.fftpack 示例具有整数个信号周期，其中日期和频率取自 FFT 理论。

编码：

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

# 1. Linspace
N = 600
# Sample spacing
tmax = 3/4
T = tmax / N # =1.0 / 800.0
x1 = np.linspace(0.0, N*T, N)
y1 = np.sin(50.0 * 2.0*np.pi*x1) + 0.5*np.sin(80.0 * 2.0*np.pi*x1)
yf1 = scipy.fftpack.fft(y1)
xf1 = np.linspace(0.0, 1.0/(2.0*T), N//2)

# 2. Integer number of periods
tmax = 1
T = tmax / N # Sample spacing
x2 = np.linspace(0.0, N*T, N)
y2 = np.sin(50.0 * 2.0*np.pi*x2) + 0.5*np.sin(80.0 * 2.0*np.pi*x2)
yf2 = scipy.fftpack.fft(y2)
xf2 = np.linspace(0.0, 1.0/(2.0*T), N//2)

# 3. Correct positioning of dates relatively to FFT theory ('arange' instead of 'linspace')
tmax = 1
T = tmax / N # Sample spacing
x3 = T * np.arange(N)
y3 = np.sin(50.0 * 2.0*np.pi*x3) + 0.5*np.sin(80.0 * 2.0*np.pi*x3)
yf3 = scipy.fftpack.fft(y3)
xf3 = 1/(N*T) * np.arange(N)[:N//2]

fig, ax = plt.subplots()
# Plotting only the left part of the spectrum to not show aliasing
ax.plot(xf1, 2.0/N * np.abs(yf1[:N//2]), label='fftpack tutorial')
ax.plot(xf2, 2.0/N * np.abs(yf2[:N//2]), label='Integer number of periods')
ax.plot(xf3, 2.0/N * np.abs(yf3[:N//2]), label='Correct positioning of dates')
plt.legend()
plt.grid()
plt.show()

输出：

就像这里一样，即使使用整数个周期，一些扩散仍然存在。这种行为是由于 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。

我只是得出结论，使用示例应该替换为以下代码（在我看来，这不那么具有误导性）：

import numpy as np
from scipy.fftpack import fft

# Number of sample points
N = 600
T = 1.0 / 800.0
x = T*np.arange(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 = 1/(N*T)*np.arange(N//2)
import matplotlib.pyplot as plt
plt.plot(xf, 2.0/N * np.abs(yf[0:N//2]))
plt.grid()
plt.show()

输出（第二个尖峰不再扩散）：

我认为这个答案仍然带来了一些关于如何正确应用离散傅立叶变换的额外解释。显然，我的回答太长了，而且总是有其他的事情要说（例如，ewerlopes 简要地谈到了混叠，关于窗口化可以说很多），所以我会停下来。

我认为在应用离散傅立叶变换时深入理解它的原理非常重要，因为我们都知道很多人在应用它时会在这里和那里添加因子以获得他们想要的东西。