scipy 在Python中绘制快速傅立叶变换

aiazj4mn  于 2022-11-10  发布在  Python
关注(0)|答案(7)|浏览(194)

我可以访问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/0WhjjMkbhttp://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]

mrphzbgm

mrphzbgm1#

因此,我在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)

siotufzp

siotufzp2#

关于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()

这应该能解决你的问题。

t30tvxxf

t30tvxxf3#

高尖峰是由于信号的DC(不变,即freq = 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]))

将显示:

am46iovg

am46iovg4#

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

sauutmhj

sauutmhj5#

作为对上述答案的补充,我想指出的是,通常情况下,调整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')

输出曲线:

xqnpmsa8

xqnpmsa86#

我写这个额外的答案是为了解释使用FFT时尖峰扩散的起源,特别是讨论我在某些地方不同意的scipy.fftpack教程。
在本例中,记录时间为tmax=N*T=0.75。信号为sin(50*2*pi*x) + 0.5*sin(80*2*pi*x)。频率信号应包含频率为5080的两个尖峰,振幅为10.5。但是,如果分析的信号不具有整数个周期,则由于信号截断,可能出现扩散:

  • 峰值1:50*tmax=37.5 =〉频率50不是1/tmax =〉的倍数,由于在该频率处的信号截断,存在扩散
  • 峰值2:80*tmax=60 =〉频率801/tmax =〉的倍数,由于在该频率处的信号截断,因此没有扩散

以下代码分析的信号与教程(sin(50*2*pi*x) + 0.5*sin(80*2*pi*x))中的信号相同,但略有不同:
1.原始的scipy.fftpack示例。
1.原始scipy.fftpack示例具有整数个信号周期(tmax=1.0代替0.75以避免截断扩散)。
1.原始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形式还导致尖峰的偏移,所述尖峰位于比它们应该在的频率稍高的频率处,如在第一幅图中可以看到的,其中尖峰在频率为5080
我将得出结论,应该用下面的代码替换用法示例(在我看来,这不容易引起误解):

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简要地谈到了aliasing,还有很多关于窗口的内容),所以我就不多说了。
我认为在应用离散傅里叶变换时,深入理解它的原理是非常重要的,因为我们都知道,很多人在应用它时,为了得到他们想要的东西,到处添加因子。

mrwjdhj3

mrwjdhj37#

本页中已经有了很好的解决方案,但所有的解决方案都假设数据集是均匀/均匀采样/分布的。我将尝试提供一个更一般的随机采样数据的示例。我还将使用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笔记本here

相关问题