使用NumPy中的FFT时遇到的频率分辨率问题

Question

我用Tektronix示波器进行信号采集，获得了一万个测量点（几个信号周期），需要对这组数据进行频率分析。我的信号为8MHz正弦波。当我使用SciPy或NumPy时，得到的结果相同 - 频率分布太宽。两个值之间的距离为500kHz，最高频率为2.5GHz（荒谬）。当我想要测量8MHz左右的频带时，只能得到7.5、8.0和8.5 MHz的确切值。我试图改变由(x[1]-x[0])确定的样本间隔，但没有得到更好的结果。

def CalculateFFT(t_val,p_val):
    x = t_val #Two parameters: [x,y] values
    y = lambda x: p_val
    com_signal = y(x) # Combined signal
    FFT_val = abs(scipy.fft(com_signal))
    freq_val = scipy.fftpack.fftfreq(len(com_signal), x[1]-x[0])
    spec_val = 20*scipy.log10(FFT_val)
    return freq_val, spec_val

- akson128

---

你的测量周期应该比几个信号周期长得多，以获得更准确的频率分辨率。请问测量的采样频率是多少？ - Jan Kuiken

不要误会，你确定你完全理解如何设计DFFT的最佳输入以及如何解释DFFT的输出（这两者都不是轻松的）吗？也许你想阅读一下我的论文中关于FFT的部分：http://gehrcke.de/files/stud/gehrcke_MScThesis_magnetic_particle_imaging.pdf（突然间我也不知道问题出在哪里）。 - Dr. Jan-Philip Gehrcke

谢谢Jan-Philip Gehrcke，这对我很有帮助（而且也是一个不错的论文话题）。我进行了一些额外的模拟，并成功地发现，如果我在恒定数据集（10k个数据点）中设置更多的信号周期，通过改变时间窗口大小，我能得到更准确的频率值。 - akson128

Answer 1

值得深入阅读DFFT（离散傅里叶变换）的工作原理，但您应该始终记住以下公式。对于具有n个点和最大时间Tmax的时间序列，时间分辨率由dt = Tmax / n给出。

DFFT将产生n个点，其中

Fmax = 1 / dt

dF = 1 / Tmax

您似乎暗示最大频率足够（因此时间分辨率可以），但是频率分辨率不够好：您需要以相同的时间分辨率收集更多数据。

- Greg

Answer 2

如果(1)采样时间太短，(2)需要更高的估计频率精度，并且(3)您知道您的信号是正弦波，那么您可以将该信号拟合成正弦波。就像在如何使用pylab和numpy将正弦曲线拟合到我的数据？中所述一样，唯一的例外是需要添加频率。 这里是一个频率约为8 MHz 的示例图：

以下是示例代码：

""" Modified from https://dev59.com/XmQn5IYBdhLWcg3wmIAs#16716964 """
from numpy import sin, linspace, pi,average;
from pylab import plot, show, title, xlabel, ylabel, subplot, scatter
from scipy import fft, arange, ifft
import scipy
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import leastsq

ff = 8e6;   # frequency of the signal
Fs = ff*128;  # sampling rate
Ts = 1.0/Fs; # sampling interval

t = arange(0,((1/ff)/128)*(128)*5,Ts) # time vector
A = 2.5;

ff_0 = 8.1456e6
y = A*np.sin(2*np.pi*ff_0*t+15.38654*pi/180) + np.random.randn(len(t))/5

guess_b = 0
guess_a = y.std()*2**0.5;
guess_c = 10*pi/180
guess_d = ff*0.98*2*pi

fig = plt.figure(facecolor="white")
plt.plot(t,y,'.', label='Signal Fred. %0.4f Hz'%(ff_0/1e6))
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.grid(alpha=0.5);

optimize_func = lambda x: (x[0]*np.sin(x[2]*t+x[1]) - y);
est_a,  est_c, est_d = leastsq(optimize_func, [guess_a, guess_c, guess_d])[0]
data_fit = est_a*np.sin(est_d*t+est_c) ;
plt.plot(t,data_fit,label='Fitted Est. Freq. %0.4f Hz'%(est_d/(2*pi)/1e6))
plt.legend()
plt.tight_layout();
plt.show();

fig.save("sinfit.png")

- user6036470