在 Python 中绘制快速傅里叶变换

Posted 2023-03-12

【中文标题】在 Python 中绘制快速傅里叶变换

【英文标题】：Plotting a fast Fourier transform in Python

【发布时间】：2014-11-02 07:34:48

【问题描述】：

我可以访问 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/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]。

【问题讨论】：

向我们展示您的尝试、失败的原因以及您正在使用的示例。

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

这是一个很好的例子，但问题到底出在哪里？该代码对我很有用。剧情根本就没有出现吗？

即，您使用的是哪种参数（我们至少需要查看您的一些数据）

我添加了 x 和 y 轴的 pastebin，x 数据以秒为单位，y 数据只是传感器读数。当我将这些数据列表放入 fft 示例时，它只是在零处出现了巨大的峰值

【参考方案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 的峰值正是我所期望的。那么有什么问题呢？

为了响应发布的原始数据和 cmets

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

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)

【讨论】：

不是示例错误，而是我不知道如何将其应用到我的数据中。

@user3123955，对。这就是为什么我们需要查看您的数据以及如果我们要帮助您，它是如何失败的。

@user3123955 那么您希望任何 FFT 算法对此做什么？你需要清理你的数据。

语句 [:N/2] 应为 [:N//2] 以避免出现弃用警告。浮点值不应该是索引。在 2.x N/2 中创建了一个 int。在 3.x N/2 中创建一个浮点数。 N//2 创建 2.x 中预期的 int

@PaulH 频率50 Hz 的幅度不应该是1，频率80 Hz 的幅度不应该是0.5？

【参考方案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()

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

【讨论】：

我已经插值了我的数据以获得均匀的间距，你能告诉我 fftfreq 的确切作用吗？为什么它需要我的 x 轴？为什么要绘制 Y 的绝对值和角度？角度是相位吗？相位与什么有关？当我对我的数据执行此操作时，它只是在 0Hz 处有一个巨大的峰值，并且很快就结束了，但是我给它提供的数据没有恒定的偏移量（我对边缘为 0.15 Gz 到 12Hz 的数据进行了大的带通）为了摆脱恒定的偏移量，我的数据无论如何都不应大于 4 Hz，因此频段应该会让我丢失信息）。

1. fftfreq 为您提供与您的数据相对应的频率分量。如果您绘制freq，您将看到 x 轴不是一个不断增加的函数。您必须确保在 x 轴上有正确的频率分量。可以看说明书：docs.scipy.org/doc/numpy/reference/generated/…

2.大多数人都喜欢查看 fft 的幅度和相位。相位信息会告诉你什么很难用一句话解释，但我只能说，当你组合信号时它是有意义的。当您组合同相的相同频率的信号时，它们会放大，而当它们异相 180 度时，它们会衰减。当您设计放大器或任何有反馈的东西时，这一点变得很重要。

3.通常，您的最低频率实际上将具有零相位，并且参考了这一点。当信号通过您的系统时，每个频率都以不同的速度移动。这是相速度。相位图为您提供了这些信息。我不知道你正在使用什么系统，所以不能给你一个明确的答案。对于此类问题，最好阅读反馈控制、模拟电子学、数字信号处理、电磁场理论等，或者更适合您的系统的内容。

4.与其使用 your 数据，不如从生成自己的信号开始：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。您将对您的编程充满信心。然后你可以判断你的结果的真实性。如果您尝试分析的信号是您使用过的第一个信号，那么您总是会觉得自己做错了什么......

【参考方案3】：

您的高尖峰是由于信号的 DC（不变，即频率 = 0）部分造成的。这是一个规模问题。如果您想查看非 DC 频率内容，为了可视化，您可能需要从信号 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]))

将显示：

【讨论】：

是的，单位是赫兹。在代码中，xf 的定义将 fft bin 映射到频率。

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

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

【参考方案4】：

作为对已经给出的答案的补充，我想指出，通常对于 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')

输出图：

【参考方案5】：

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

【讨论】：

这看起来非常接近我对音乐频段显示的需求：每 33 毫秒（每秒 30 帧）拍摄一次系统声音快照。将频率分为 3、5、7、9 或 11 个频段。计算 100% 的频带幅度百分比。不需要matplotlib，因为我将在 tkinter 画布上绘制 LED。您能指出这样的答案吗，或者如果我问一个新问题，您会好心回答吗？谢谢，

【参考方案6】：

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

派克 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))，但略有不同：

原始 scipy.fftpack 示例。 带有整数个信号周期的原始 scipy.fftpack 示例（tmax=1.0 代替 0.75 以避免截断扩散）。 原始 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 talked briefly 关于aliasing 可以说很多关于windowing），所以我会停下来。

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

【参考方案7】：

此页面上已经有很好的解决方案，但都假设数据集是均匀/均匀采样/分布的。我将尝试提供一个更一般的随机抽样数据示例。我也会以this MATLAB tutorial为例：

添加所需的模块：

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)

P.S.我终于有时间实现一个更规范的算法来获得不均匀分布数据的傅里叶变换。您可能会看到代码、描述和示例 Jupyter notebook here。

【讨论】：

我在文档中没有看到任何建议 resample 处理非均匀采样时间的内容。它确实接受时间参数（示例中未使用），但这似乎也假设了统一的采样时间。

@user2699 this example 可能会有所帮助

'scipy.signal.resample` 使用 FFT 方法对数据进行重新采样。用它对非均匀数据重新采样以获得均匀的 FFT 是没有意义的。

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

您给出的所有方法都有优点和缺点（尽管请注意sklearn.utils.resample 不执行插值）。如果您想讨论可用于查找不规则采样信号的频率的选项，或者不同类型插值的优点，请开始另一个问题。两者都是有趣的主题，但远远超出了如何绘制 FFT 的答案范围。