基于数据集大小的 FFT 偏移频率？

Question

我正在努力从给定的数据集中查找频率，我正在努力理解 np.fft.fft() 的工作原理。我以为我有一个工作脚本，但遇到了一个我无法理解的奇怪问题。

我有一个大致为正弦曲线的数据集，我想了解信号由哪些频率组成。一旦我进行了 FFT，我就得到了这个情节：

但是，当我采用相同的数据集，将其切成两半并绘制相同的东西时，我得到了：

我不明白为什么频率从 144kHz 下降到 128kHz，这在技术上应该是相同的数据集，但长度更小。

我可以确认几件事：

1. 数据点之间的步长 0.001
2. 我尝试过插值，但运气不佳。
3. 如果我对数据集的后半部分进行切片，我也会得到不同的频率。
4. 如果我的数据集确实由 128 和 144kHz 组成，那么为什么第一个图中没有出现 128 峰值？

更令人困惑的是，我正在运行一个没有问题的纯正弦波脚本：

T = 0.001
fs = 1 / T
def find_nearest_ind(data, value):
    return (np.abs(data - value)).argmin()
x = np.arange(0, 30, T)
ff = 0.2
y = np.sin(2 * ff * np.pi * x)
x = x[:len(x) // 2]
y = y[:len(y) // 2]

n = len(y) # length of the signal
k = np.arange(n)
T = n / fs
frq = k / T * 1e6 / 1000 # two sides frequency range
frq = frq[:len(frq) // 2] # one side frequency range

Y = np.fft.fft(y) / n # dft and normalization
Y = Y[:n // 2]

frq = frq[:50]
Y = Y[:50]

fig, (ax1, ax2) = plt.subplots(2)
ax1.plot(x, y)
ax1.set_xlabel("Time (us)")
ax1.set_ylabel("Electric Field (V / mm)")
peak_ind = find_nearest_ind(abs(Y), np.max(abs(Y)))
ax2.plot(frq, abs(Y))
ax2.axvline(frq[peak_ind], color = 'black', linestyle = '--', label = F"Frequency = {round(frq[peak_ind], 3)}kHz")
plt.legend()
plt.xlabel('Freq(kHz)')
ax1.title.set_text('dV/dX vs. Time')
ax2.title.set_text('Frequencies')
fig.tight_layout()
plt.show()

Answer 1

这里是您的代码的细分，以及一些改进建议和额外的解释。仔细研究它会告诉你发生了什么。你得到的结果完全是预期的。最后我会提出一个通用的解决方案。

首先正确设置您的单位。我假设您正在处理秒，而不是微秒。只要保持一致，您就可以稍后进行调整。

确定采样的周期和频率。这意味着 FFT 的奈奎斯特频率将为 500Hz：

T = 0.001      # 1ms sampling period
fs = 1 / T     # 1kHz sampling frequency

制作一个 30e3 个点的时域。半域将包含 15000 个点。这意味着频率分辨率为 500Hz / 15k = 0.03333Hz。

x = np.arange(0, 30, T)   # time domain
n = x.size                # number of points: 30000

在做任何事情之前，我们可以在这里定义我们的时域。我更喜欢比您使用的更直观的方法。这样您就不必重新定义T 或引入辅助变量k。但只要结果相同，就无所谓了：

F = np.linspace(0, 1 - 1/n, n) / T    # Notice F[1] = 0.03333, as predicted

现在定义信号。你选择了ff = 0.2。注意 0.2Hz。 0.2 / 0.03333 = 6，因此您会期望在 bin 索引 6 (F[6] == 0.2) 中看到峰值。为了更好地说明正在发生的事情，让我们以ff = 0.22 为例。这会将频谱泄漏到相邻的 bin 中。

ff = 0.22
y = np.sin(2 * np.pi * ff * x)

现在进行 FFT：

Y = np.fft.fft(y) / n
maxbin = np.abs(Y).argmax()  # 7
maxF = F[maxbin]             # 0.23333333: This is the nearest bin

由于您的频率区间为 0.03Hz，因此您可以预期的最佳分辨率为 0.015Hz。对于分辨率低得多的真实数据，误差要大得多。

现在让我们看看将数据大小减半时会发生什么。其中，频率分辨率变得更小。现在你有一个 500Hz 的最大频率分布在 7.5k 样本上，而不是 15k：分辨率下降到每个 bin 0.066666Hz：

n2 = n // 2                               # 15000
F2 = np.linspace(0, 1 - 1 / n2, n2) / T   # F[1] = 0.06666
Y2 = np.fft.fft(y[:n2]) / n2

看看频率估计会发生什么：

maxbin2 = np.abs(Y2).argmax()  # 3
maxF2 = F2[maxbin2]            # 0.2: This is the nearest bin

希望您能看到这如何应用于您的原始数据。在完整的 FFT 中，完整数据的每个 bin 的分辨率约为 16.1，半数据的分辨率约为 32.2kHz。因此，您的原始结果在右峰值的 ~±8kHz 范围内，而第二个结果在 ~±16kHz 范围内。因此，真实频率在 136kHz 和 144kHz 之间。另一种看待它的方法是比较你给我看的垃圾箱：

满：128.7 144.8 160.9 一半：96.6 128.7 160.9

当您取出恰好一半的数据时，您会删除每隔一个频率区间。如果您的峰值最初最接近 144.8kHz，而您放弃了那个 bin，它将最终变为 128.7 或 160.9。

注意：根据您显示的 bin 编号，我怀疑您对 frq 的计算有一点偏差。注意我的linspace 表达式中的1 - 1/n。你需要它来获得正确的频率轴：最后一个 bin 是 (1 - 1/n) / T，而不是 1 / T，不管你如何计算它。

那么如何解决这个问题呢？最简单的解决方案是对峰值周围的三个点进行抛物线拟合。当您正在寻找基本完美的正弦曲线时，这通常是对数据中真实频率的足够好的估计。

def peakF(F, Y):
    index = np.abs(Y).argmax()
    # Compute offset on normalized domain [-1, 0, 1], not F[index-1:index+2]
    y = np.abs(Y[index - 1:index + 2])
    # This is the offset from zero, which is the scaled offset from F[index]
    vertex = (y[0] - y[2]) / (0.5 * (y[0] + y[2]) - y[1])
    # F[1] is the bin resolution
    return F[index] + vertex * F[1]

如果您想知道我是如何得到抛物线公式的：我用x = [-1, 0, 1] 和y = Y[index - 1:index + 2] 求解了系统。矩阵方程为

[(-1)^2 -1  1]   [a]   Y[index - 1]
[   0^2  0  1] @ [b] = Y[index]
[   1^2  1  1]   [c]   Y[index + 1]

使用归一化域计算偏移量并在之后进行缩放几乎总是比使用 F[index - 1:index + 2] 中的任何巨大数字更稳定。

您可以将示例中的结果插入到这个函数中，看看它是否有效：

>>> peakF(F, Y)
0.2261613409657391
>>> peakF(F2, Y2)
0.20401580936430794

如您所见，抛物线拟合带来了改善，但幅度很小。但是，通过更多样本来提高频率分辨率是无法替代的！