在 numpy 中使用 FFT 的频率分辨率问题答案

【问题标题】：Frequency resolution issue using FFT in numpy在 numpy 中使用 FFT 的频率分辨率问题
【发布时间】：2013-07-13 15:07:27
【问题描述】：

我使用泰克示波器来执行一些信号采集。我得到 10.000 个测量点（几个信号周期），我必须对那组数据进行频率分析。我的信号是 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

【问题讨论】：

您的测量周期应该比几个信号周期长得多以获得更准确的频率区间，测量的采样频率是多少？
无意冒犯，您是否确保您完全了解必须如何设计 DFFT 的最佳输入以及必须如何解释 DFFT 的输出（两者都不是微不足道的）？也许你可能想在这里阅读我论文的 FFT 部分：gehrcke.de/files/stud/…（我自然也不知道问题出在哪里）
谢谢 Jan-Philip Gehrcke，它很有帮助（以及不错的论文主题）。我做了一些额外的模拟，并且我设法看到通过更改时间窗口在恒定数据集（10k 点）中设置的信号周期越多，我得到的频率值就越准确。

标签： python numpy scipy frequency-analysis

【解决方案1】：

值得深入阅读 DFFT 的工作原理，但您应该始终牢记以下公式。对于具有 n 个点和最大时间 Tmax 的时间序列，时间分辨率由 dt = Tmax / n

给出

一个DFFT将产生n个点

Fmax = 1 / dt

dF = 1 / Tmax

您似乎建议最大频率足够（所以时间分辨率还可以），但频率分辨率不够好：您需要收集更多数据，同时分辨率。

【讨论】：

【解决方案2】：

如果 (1) 采样时间太短，(2) 您需要更高的估计频率精度，并且 (3) 您知道您的信号是正弦波，那么您可以将信号拟合为正弦波。比如How do I fit a sine curve to my data with pylab and numpy?，除了需要添加频率。

这是一个频率约为 8 MHz 的示例图：

下面是示例代码：

""" Modified from https://stackoverflow.com/a/16716964/6036470 """
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")

【讨论】：