在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]。

- user3123955

---

展示你尝试过的内容，它们失败的原因，以及你正在参考的示例。 - Paul H

我发布了我尝试过的示例以及我的想法，我认为我只是对如何正确绘制输出感到困惑。 - user3123955

这是一个很好的例子，但问题到底是什么？对我来说，那段代码运行得很好。难道是图表没有显示出来吗？ - Paul H

请提供您使用的参数（我们需要查看至少一些数据）。 - Paul H

我已经添加了X轴和Y轴的Pastebin，其中X数据以秒为单位，而Y数据只是传感器读数。当我将这些数据列表放入FFT示例中时，它只在零处产生了一个巨大的峰值。 - user3123955

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)

- Paul H

---

并不是这个例子有问题，而是我不知道如何将它应用到我的数据上。 - user3123955

@user3123955，没错。这就是为什么如果你想要帮助，我们需要看到你的数据以及它失败的情况。 - Paul H

@user3123955，那你希望任何FFT算法对此做些什么呢？你需要清理你的数据。 - Paul H

语句[:N/2]应改为[:N//2]以避免弃用警告。浮点值不应作为索引。在2.x中，N/2创建了一个int。在3.x中，N/2创建了一个float。N//2创建了2.x中预期的int。 - KeithSmith

@PaulH，频率为50 Hz的振幅应该是1，频率为80 Hz的振幅应该是0.5吗？ - Furqan Hashim

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]）。在这种情况下，你可以直接使用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()

这应该解决你的问题。

- ssm

---

我已经对我的数据进行了插值以获得均匀间距，你能告诉我fft频率到底是做什么用的吗？为什么需要我的x轴？你为什么要绘制Y的绝对值和角度？角度是否表示相位？相位与什么有关？当我使用我的数据进行这样的操作时，它只在0Hz处有一个巨大的峰值，然后很快就衰减了，但是我正在输入没有常量偏移的数据（我对数据进行了较大的带通滤波，边缘为0.15 Gz到12 Hz，以消除恒定的偏移，我的数据不应大于4 Hz，因此该带通滤波应使我丢失信息）。 - user3123955

fftfreq函数能够给出与你的数据相对应的频率分量。如果绘制freq，你会发现x轴不是一个持续增加的函数。你需要确保x轴上有正确的频率分量。你可以查看手册：http://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.fftfreq.html - ssm

大多数人喜欢查看FFT的幅度和相位。很难用一句话解释相位信息会告诉您什么，但我能说的是，当您合并信号时，它是有意义的。当您合并具有相同频率且处于同相位的信号时，它们会放大，而当它们相位相差180度时，它们将减弱。这在设计放大器或任何具有反馈的设备时非常重要。 - ssm

通常，您的最低频率将几乎没有相位，并且它是以此为参考的。当信号通过您的系统移动时，每个频率都以不同的速度移动，这就是相速度。相图提供了这些信息。我不知道您使用的系统是什么，因此无法给出确定的答案。对于此类问题，最好阅读有关反馈控制、模拟电子学、数字信号处理、电磁场理论等方面的资料，或者更具体地针对您的系统。 - ssm

为了避免使用您的数据，请从生成您自己的信号开始：t = linspace(0, 10, 1000); ys = [ (1.0/i)*sin(i*t) for i in arange(10)]; y = reduce(lambda m, n: m+n, ys)。然后绘制每个ys和总y，并获取每个分量的FFT。这样可以增强您对编程的信心。接下来，您可以判断结果的真实性。如果您尝试分析的信号是您第一次进行FFT分析的信号，那么您始终会感到自己在做错事... - ssm

好的答案。给出振幅和相位谱。 - Paulo Carvalho

Answer 3

我已经编写了一个处理绘制实信号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)
        + 10 * np.sin(2 * np.pi * (f0 * 2 + 0.5) * t)  # <--- not sampled on grid so the peak will not be actual height
    )
    # Result in frequencies
    fftPlot(sig, dt=dt)
    # Result in samples (if the frequencies axis is unknown)
    fftPlot(sig)

- YoniChechik

---

这看起来非常接近我需要的音乐频率带显示：每33毫秒（每秒30帧）拍摄系统声音快照。将频率分成3、5、7、9或11个频带。计算100%的带幅度百分比。不需要使用matplotlib，因为我会在tkinter画布上绘制LED灯。您能否指向一个类似这样的答案，或者如果我提出一个新问题，您是否可以友好地回答？谢谢。 - WinEunuuchs2Unix

@YoniChechik：感谢分享！我想知道这个能否用于实时信号，不是等间隔的？ - shoggananna

@user9106985 是的，需要进行必要的修改。尝试阅读有关STFT的内容。 - YoniChechik

Answer 4

你的高峰是由信号的直流（非变化，即freq = 0）部分引起的。这是一个规模问题。如果你想看到非直流频率内容以进行可视化，你可能需要从信号FFT的偏移1而不是偏移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]))

将显示：

- hesham_EE

---

是的，它是以赫兹为单位的。在代码中，“xf”的定义将FFT bin映射到频率。 - hesham_EE

不错！那么y轴呢？振幅呢？非常感谢，hesham_EE。 - vicemagui

是的，y轴是复fft的绝对值。请注意使用np.abs()函数。 - hesham_EE

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

输出图形：

- ewerlopes

Answer 6

我写这个额外的答案来解释在使用FFT时扩散峰的起源，特别是讨论我在某些地方不同意的scipy.fftpack教程。

在这个例子中，记录时间为tmax=N*T=0.75。信号为sin(50*2*pi*x) + 0.5*sin(80*2*pi*x)。频率信号应该包含两个幅值分别为1和0.5，频率分别为50和80的峰。然而，如果被分析的信号没有一个整数周期，由于信号被截断，就会出现扩散：

• 峰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. 源自FFT理论的具有整数信号周期的原始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是采样周期，信号的总持续时间为。请注意，我们停止于。
• 相关频率为，其中是采样频率。所有信号的谐波应为采样频率的倍数，以避免扩散。

在上面的示例中，可以看到使用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简要介绍了混叠，还有许多关于窗函数的内容可以说），所以我就停止了。

我认为，在应用离散傅里叶变换时深入理解其原理非常重要，因为我们都知道很多人在应用它时会添加各种因素，以获得他们想要的结果。

- bousof

Answer 7

这个页面上已经有很好的解决方案，但是所有的解决方案都假设数据集是均匀/平均采样/分布的。我将尝试提供更一般的随机抽样数据的示例。我还将使用这个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")

Now calculating the 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. 我终于有时间实现了一种更加规范的算法来获取不均匀分布数据的傅里叶变换。您可以在这里查看代码、描述和示例Jupyter笔记本。

- Foad S. Farimani

---

scipy.signal.resample 使用FFT方法对数据进行重采样。将非均匀数据重采样为均匀FFT是没有意义的，因此不应使用该函数进行操作。 - user2699

@user2699 scipy.signal.resample 只是其中一个可用的选项，不一定是最好的。我会检查其他选项，进行比较，并编辑我的帖子。感谢您的评论。 - Foad S. Farimani

@user2699 看来我在这里太天真了。已经有一些可用的库了：1. nfft 库，它似乎是 NFFT 的包装器；2. pyNFFT 和 3. PyNUFFT。 - Foad S. Farimani

scipy.interpolate，sklearn.utils.resample，pandas.DataFrame.resample，numpy.interp 是其他可用的插值/重采样选项。 - Foad S. Farimani

你提供的所有方法都有优点和缺点（请注意，sklearn.utils.resample 不执行插值）。如果您想讨论查找不规则采样信号频率的选项或不同类型插值的优点，请提出另一个问题。这两个主题都很有趣，但超出了如何绘制FFT的答案范围。 - user2699