Scipy/Numpy FFT 频率分析

Question

Scipy/Numpy FFT 频率分析

33

我正在寻找如何将通过scipy.fftpack.fftfreq获取的fft的频率轴转换为赫兹频率，而不是二进制或分数二进制。我尝试了以下代码来测试FFT：

t = scipy.linspace(0,120,4000)
acc = lambda t: 10*scipy.sin(2*pi*2.0*t) + 5*scipy.sin(2*pi*8.0*t) + 2*scipy.random.random(len(t))

signal = acc(t)

FFT = abs(scipy.fft(signal))
FFT = scipy.fftpack.fftshift(FFT)
freqs = scipy.fftpack.fftfreq(signal.size)

pylab.plot(freqs,FFT,'x')
pylab.show()

采样率应为4000个样本/120秒=33.34个样本/秒。

信号具有2.0赫兹的信号、8.0赫兹的信号和一些随机噪声。

我进行FFT处理，获取频率并绘制图形。这些数字相当荒谬。如果我将频率乘以33.34（采样频率），那么我会在大约8赫兹和15赫兹处得到峰值，这似乎是错误的（此外，频率应该相差4倍而不是2倍！）。

对此您有什么想法吗？

- nathan lachenmyer

4个回答

11

scipy.fftpack.fftfreq(n, d) 可以直接得到频率。如果您将d=1/33.34，这将告诉您每个fft点的频率（以Hz为单位）。

- tom10

6

每个频率段的宽度为（采样频率/频率段数）。更加根本的问题是，你的采样率不足以捕获感兴趣的信号。你的采样率为8.3赫兹；你需要至少16赫兹才能捕获8赫兹的输入音调。注1：对于所有数字信号处理专家来说，实际上与其相关的是带宽而不是最大频率。但我假设原作者不想进行欠采样数据采集。

- Oliver Charlesworth

我使用了4000个样本，持续120秒--这不是33.3赫兹吗？这应该足够了，但数字仍然不准确... - nathan lachenmyer

@asymptoticdesign：啊，好的。你最初的问题说的是1000。是的，那应该足够了。这个能量出现在哪个箱子索引中？ - Oliver Charlesworth

-2

你的方程式出了问题。

fs = 33.33
df1 = 2*pi * (2.0/fs)
df2 = 2*pi * (5.0/fs)
x = [10*sin(n*df1) + 5*sin(n*df2) + 2*random.random() for n in range(4000)]

这将为您提供4000个样本，采样率为33.33 Hz，代表120秒的数据。

现在进行FFT。Bin 0将保存DC结果。Bin 1将是33.33，bin 2将是66.66，以此类推。

编辑：我忘记提到，由于您的采样率为33.33 Hz，因此可以表示的最大频率将是fs/2，即16.665 Hz。

- sizzzzlerz

2

-1：不是的。总带宽为33.33Hz，而不是每个频段的宽度。 - Oliver Charlesworth

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- HYRY · Accepted Answer

我认为你不需要使用fftshift()函数，而是可以直接将采样周期传递给fftfreq()函数：

import scipy
import scipy.fftpack
import pylab
from scipy import pi
t = scipy.linspace(0,120,4000)
acc = lambda t: 10*scipy.sin(2*pi*2.0*t) + 5*scipy.sin(2*pi*8.0*t) + 2*scipy.random.random(len(t))

signal = acc(t)

FFT = abs(scipy.fft(signal))
freqs = scipy.fftpack.fftfreq(signal.size, t[1]-t[0])

pylab.subplot(211)
pylab.plot(t, signal)
pylab.subplot(212)
pylab.plot(freqs,20*scipy.log10(FFT),'x')
pylab.show()

从图中可以看出，有两个峰值分别在2Hz和8Hz。

enter image description here