Scipy/Numpy FFT频率分析

Question

我正在寻找如何将fft中的频率轴(通过scipy.fftpack.fftfreq获取)转换为赫兹频率,而不是箱子或分数箱.

我尝试在下面编写代码来测试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 Hz信号,8.0 Hz信号和一些随机噪声.

我采用FFT,抓取频率并绘制它.这些数字非常荒谬.如果我将频率乘以33.34(采样频率),那么我会得到大约8 Hz和15 Hz的峰值,这似乎是错误的(同样,频率应该是4倍,而不是2!).

对我在这里做错了什么的想法？

Answer 1

我认为你不需要做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有两个峰值.

Answer 2

scipy.fftpack.fftfreq(n,d)直接给你频率.如果设置d=1/33.34,这将告诉您fft每个点的频率(Hz).

Answer 3

每个bin的频率宽度是(sampling_freq/num_bins).

更基本的问题是您的采样率不足以满足您感兴趣的信号.您的采样率为8.3 Hz; 您需要至少16Hz才能捕获8Hz输入音.¹

_{1.致所有DSP专家; 我知道实际上BW是相关的,而不是最大频率.但我假设OP不想做欠采样数据采集.}