如何在快速傅立叶变换中正确缩放频率轴? [英] How to properly scale frequency axis in Fast Fourier Transform?

问题描述

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

I am trying some sample code taking the FFT of a simple sinusoidal function. Below is the code

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处有两个尖峰.但是，当我运行此代码时，得到以下图.

From the function that is given, , it is clear there should be two spikes at frequencies 1 and 5. However, when I run this code, I get the following plot.

                                   

                                         

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

Clearly the spikes are not where they should be. Additionally, I have noticed that the frequency scaling is sensitive to the number of points N as well as the interval limits that I make limit. As an example, setting N = 2048 gives the following plot.

                                   

                                         

如您所见，尖峰的位置已更改.现在，保持N = 1024并将其设置为limit = 100也会更改结果.

As you can see, the locations of the spikes have changed. Now keeping N = 1024 and setting limit = 100 also changes the result.

                                   

                                         

我该如何做才能使频率轴始终正确缩放?

推荐答案

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 执行实际的(快速)傅立叶变换.它对输入采样进行相同的假设，即等距，并以与fftfreq相同的顺序输出傅立叶分量.它不在乎实际的频率值:采样间隔不会作为参数传递.

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.

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

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.

有一个突出的频率:f = 0.它是信号平均值(从零开始的偏移量)的量度.在fft返回的频谱和fftfreq的频率范围中，它位于第一个数组索引处. 如果绘制两个半部，则可能需要左右移动频谱，以使负半部位于零分量的左侧，正半部位于零分量的右侧，这意味着所有值都在升序并准备绘制.

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 确实可以那.但是，如果您只绘制了频谱的一半，那么您可能根本不会理会它.尽管您使用 if 进行操作，但必须同时移动两个 数组:频率和傅立叶分量.在您的代码中，您仅移动了频率.这就是峰值最终出现在频谱的错误一侧的原因:您将傅立叶分量绘制成相对于频率的正一半相对于负一半，因此右侧的峰值实际上意味着接近于零，而不是在远端.

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.

您实际上并不需要依赖那些在频率上运行的功能.仅根据fftfreq的文档即可生成其范围:

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

请注意，Numpy还提供了专用功能， rfft rfftfreq ，用于实值数据的常见用例.

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

这篇关于如何在快速傅立叶变换中正确缩放频率轴?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！