如何在快速傅立叶变换中正确缩放频率轴？

Question

我正在尝试一些采用简单正弦函数的 FFT 的示例代码。下面是代码

import numpy as np
from matplotlib import pyplot as plt

N = 1024
limit = 10
x = np.linspace(-limit, limit, N)
dx = x[1] - x[0]
y = np.sin(2 * np.pi * 5 * x) + np.sin(2 * np.pi * x)
Y = np.abs(np.fft.fft(y) ** 2)
z = fft.fftshift(np.fft.fftfreq(N, dx))
plt.plot(z[int(N/2):], Y[int(N/2):])
plt.show()

根据给定的函数 $y = \sin(10\pi x) + \sin(2\pi x)$ ，很明显在频率 1 和 5 处应该有两个尖峰。但是，当我运行这段代码时，我得到了以下图表。

很明显尖峰不在它们应该在的位置。此外，我注意到频率缩放对点数 N 以及我设置的间隔限制 limit 很敏感。例如，设置 N = 2048 给出以下图。

如您所见，尖峰的位置发生了变化。现在保持 N = 1024 并设置 limit = 100 也会更改结果。

我怎样才能使频率轴始终正确缩放？

Answer 1

fftfreq returns 频率范围按以下顺序：正频率从最低到最高，然后负频率按绝对值倒序排列。（您通常只想绘制一半，就像您在代码中所做的那样。）请注意函数实际上需要知道的数据很少：仅需要样本数量及其在时域中的间距。

fft 执行实际的（快速）傅里叶变换。它对输入采样做出相同的假设，即它是等距的，并以与 fftfreq 相同的顺序输出傅里叶分量。它不关心实际频率值：采样间隔不作为参数传入。

但是它确实接受复数作为输入。在实践中，这种情况很少见。输入通常是实数样本，如上例所示。在这种情况下，傅里叶变换有一个特殊的属性：它在频域中是对称的，即 f 和 −f 具有相同的值。出于这个原因，绘制频谱的两半通常没有意义，因为它们包含相同的信息。

有一个频率很突出：f = 0。它是信号平均值的度量，它与零的偏移量。在 fft 返回的频谱和 fftfreq 返回的频率范围中，它位于第一个数组索引处。如果绘制两半，移动频谱可能是有意义的，这样负半在 zero-component 的左边，正半在它的右边，这意味着所有值都按升序排列并准备好绘制。

fftshift 正是这样做的。但是，如果您只绘制了一半的频谱，那么您可能根本就不用这样做。虽然 if 你这样做了，但你必须移动 both 数组：频率和傅里叶分量。在您的代码中，您只改变了频率。这就是峰值最终出现在频谱错误一侧的原因：您绘制了傅里叶分量，指的是频率的正半部分与负半部分，因此右侧的峰值实际上意味着接近于零，而不是在远端。

您真的不需要依赖任何在频率上运行的函数。仅根据 fftfreq 的文档生成它们的范围很简单：

from numpy.fft import fft
from numpy import arange, linspace, sin, pi as π
from matplotlib import pyplot

def FFT(t, y):
    n = len(t)
    Δ = (max(t) - min(t)) / (n-1)
    k = int(n/2)
    f = arange(k) / (n*Δ)
    Y = abs(fft(y))[:k]
    return (f, Y)

t = linspace(-10, +10, num=1024)
y = sin(2*π * 5*t) + sin(2*π * t)
(f, Y) = FFT(t, y)
pyplot.plot(f, Y)
pyplot.show()

请注意，NumPy 还提供专用函数 rfft and rfftfreq，用于 real-valued 数据的常见用例。

Answer 2

import numpy as np
from matplotlib import pyplot as plt

N = 1024
limit = 10

x = np.linspace(-limit, limit, N)
dx = x[1] - x[0]

y = np.sin(2 * np.pi * 5 * x) + np.sin(2 * np.pi * x)
yhat = np.fft.fft(y)

Y = np.abs(yhat)
freq = np.linspace(0.0, 1.0/(2*dx), N//2)

plt.plot(freq, Y[0:N//2]*(2/N))
plt.show()

如何在快速傅立叶变换中正确缩放频率轴？

How to properly scale frequency axis in Fast Fourier Transform?

python

numpy

fft

frequency