python - 使用 python 的二维 FFT 导致频率稍微偏移

Question

我知道有几个关于在 python 中使用快速傅立叶变换 (FFT) 方法的问题，但不幸的是，没有一个可以帮助我解决我的问题:

我想用 python 计算给定二维信号 f 的快速傅里叶变换，即 f(x,y)。 Python 的文档帮助很大，解决了 FFT 带来的一些问题，但与我期望它显示的频率相比，我最终还是得到了一个稍微偏移的频率。这是我的 Python 代码:

from scipy.fftpack import fft, fftfreq, fftshift
import matplotlib.pyplot as plt
import numpy as np
import math

fq = 3.0 # frequency of signal to be sampled
N = 100.0 # Number of sample points within interval, on which signal is considered
x = np.linspace(0, 2.0 * np.pi, N) # creating equally spaced vector from 0 to 2pi, with spacing 2pi/N
y = x
xx, yy = np.meshgrid(x, y) # create 2D meshgrid
fnc = np.sin(2 * np.pi * fq * xx) # create a signal, which is simply a sine function with frequency fq = 3.0, modulating the x(!) direction
ft = np.fft.fft2(fnc) # calculating the fft coefficients

dx = x[1] - x[0] # spacing in x (and also y) direction (real space)
sampleFrequency = 2.0 * np.pi / dx
nyquisitFrequency = sampleFrequency / 2.0

freq_x = np.fft.fftfreq(ft.shape[0], d = dx) # return the DFT sample frequencies 
freq_y = np.fft.fftfreq(ft.shape[1], d = dx)

freq_x = np.fft.fftshift(freq_x) # order sample frequencies, such that 0-th frequency is at center of spectrum 
freq_y = np.fft.fftshift(freq_y)

half = len(ft) / 2 + 1 # calculate half of spectrum length, in order to only show positive frequencies

plt.imshow(
    2 * abs(ft[:half,:half]) / half,
    aspect = 'auto',
    extent = (0, freq_x.max(), 0, freq_y.max()),
    origin = 'lower',
    interpolation = 'nearest',
)
plt.grid()
plt.colorbar()
plt.show()

我在运行它时得到的结果是:

现在您看到 x 方向的频率并不完全在 fq = 3 处，而是稍微向左偏移了一点。为什么是这样？我认为这与 FFT 是一种使用对称参数和

half = len(ft) / 2 + 1

用于在适当的位置显示频率。但我不太明白确切的问题是什么以及如何解决它。

编辑:我也尝试过使用更高的采样频率 (N = 10000.0)，这并没有解决问题，而是将频率稍微向右移动了太多。所以我很确定问题不在于采样频率。

注意:我知道泄漏效应会导致非物理振幅，但在这篇文章中，我主要对正确的频率感兴趣。

Answer 1

我发现了一些问题

你使用 2 * np.pi 两次，如果你想要一个很好的整数周期，你应该选择 linspace 或 arg 中的一个以正弦为弧度

另外 np.linspace 默认为 endpoint=True，给你一个额外的点 101 而不是 100

fq = 3.0 # frequency of signal to be sampled
N = 100 # Number of sample points within interval, on which signal is considered
x = np.linspace(0, 1, N, endpoint=False) # creating equally spaced vector from 0 to 2pi, with spacing 2pi/N
y = x
xx, yy = np.meshgrid(x, y) # create 2D meshgrid
fnc = np.sin(2 * np.pi * fq * xx) # create a signal, which is simply a sine function with frequency fq = 3.0, modulating the x(!) direction

您可以检查这些问题:

len(x)
Out[228]: 100

plt.plot(fnc[0])

现在修复 linspace 端点意味着您有偶数个 fft bin，因此您可以在 half calc 中删除 + 1

matshow() 似乎有更好的默认值，您的 extent = (0, freq_x.max(), 0, freq_y.max()), 在 imshow 似乎混淆了 fft bin 编号

from scipy.fftpack import fft, fftfreq, fftshift
import matplotlib.pyplot as plt
import numpy as np
import math

fq = 3.0  # frequency of signal to be sampled
N = 100  # Number of sample points within interval, on which signal is considered
x = np.linspace(0, 1, N, endpoint=False)  # creating equally spaced vector from 0 to 2pi, with spacing 2pi/N
y = x
xx, yy = np.meshgrid(x, y)  # create 2D meshgrid
fnc = np.sin(2 * np.pi * fq * xx)  # create a signal, which is simply a sine function with frequency fq = 3.0, modulating the x(!) direction

plt.plot(fnc[0])

ft = np.fft.fft2(fnc)  # calculating the fft coefficients

#dx = x[1] - x[0]  # spacing in x (and also y) direction (real space)
#sampleFrequency = 2.0 * np.pi / dx
#nyquisitFrequency = sampleFrequency / 2.0
#
#freq_x = np.fft.fftfreq(ft.shape[0], d=dx)  # return the DFT sample frequencies 
#freq_y = np.fft.fftfreq(ft.shape[1], d=dx)
#
#freq_x = np.fft.fftshift(freq_x)  # order sample frequencies, such that 0-th frequency is at center of spectrum 
#freq_y = np.fft.fftshift(freq_y)

half = len(ft) // 2  # calculate half of spectrum length, in order to only show positive frequencies

plt.matshow(
            2 * abs(ft[:half, :half]) / half,
            aspect='auto',
            origin='lower'
            )
plt.grid()
plt.colorbar()
plt.show()

放大图: