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

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()

根据给出的功能， ，很明显应该在频率1和5处有两个尖峰。但是，当我运行此代码时，得到以下图。

                                         

显然，峰值不在应有的位置。另外，我注意到频率缩放对点数N以及我设定的间隔限制很敏感limit。例如，设置N = 2048给出以下图表。

                                         

如您所见，尖峰的位置已更改。现在，保留N = 1024和设置limit = 100也会更改结果。

                                         

如何使频率轴始终正确缩放？

Answer 1

fftfreq returns the frequency range in the following order: the positive frequencies from lowest to highest, then the negative frequencies in reverse order of absolute value. (You usually only want to plot one half, as you do in your code.) Note how the function actually needs to know very little about the data: just the number of samples and their spacing in the time domain.

fft performs the actual (Fast) Fourier transformation. It makes the same assumption about the input sampling, that it's equidistant, and outputs the Fourier components in the same order as fftfreq. It doesn't care about the actual frequency values: the sampling interval is not passed in as a parameter.

It does however accept complex numbers as input. In practice, this is rare. Input is usually samples of real numbers, as in the above example. In that case, the Fourier transform has a special property: it's symmetric in the frequency domain, i.e. has the same value for f and ?f. For that reason, it often doesn't make sense to plot both halves of the spectrum, as they contain the same information.

There is one frequency that stands out: f = 0. It's a measure of the average value of the signal, its offset from zero. In the spectrum returned by fft and the frequency range from fftfreq, it's at the very first array index. If plotting both halves, it may make sense to shift the frequency spectrum around, so that the negative half is to the left of the zero-component, the positive half to its right, meaning all values are in ascending order and ready to be plotted.

fftshift does exactly that. However, if you plot only half of the spectrum anyway, then you may as well not bother doing this at all. Though if you do, you must shift both arrays: the frequencies and the Fourier components. In your code, you only shifted the frequencies. That's how the peaks ended up on the wrong side of the spectrum: You plotted the Fourier components referring to the positive half of the frequencies with respect to the negative half, so the peaks on the right are actually meant to be close to zero, not at the far end.

You don't really need to rely on any of those functions operating on the frequencies. It is straightforward to generate their range based on the documentation of fftfreq alone:

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()

Note that Numpy also offers dedicated functions, rfft rfftfreq, for the common use case of real-valued data.