在 FFT 的频率上计算 dFT

Question

我正在计算函数 f(x) 在 x_i、i=0,1,...,N（已知 dx）处采样的 dFT，频率为 u_j, j=0,1,... ,N 其中 u_j 是 np.fft.fftfreq(N, dx) 生成的频率，并将其与 np.fft.fft(f(x)) 的结果进行比较。我发现这两个不同意...

我错过了什么吗？根据定义，它们不应该相同吗？ （当我查看 dFT/FFT 的图像部分时，差异甚至更糟）。

我附上了我使用的脚本，它生成了比较 dFT 和 FFT 的实部和图像部分的图。

import numpy as np
import matplotlib.pyplot as plt
from astropy import units

def func_1D(x, sigma_x):
    return np.exp(-(x**2.0 / (2.0 * sigma_x**2)))


n_pixels = int(2**5.0)
pixel_scale = 0.05 # units of arcsec

x_rad = np.linspace(
    -n_pixels * pixel_scale / 2.0 * (units.arcsec).to(units.rad) + pixel_scale / 2.0 * (units.arcsec).to(units.rad),
    +n_pixels * pixel_scale / 2.0 * (units.arcsec).to(units.rad) - pixel_scale / 2.0 * (units.arcsec).to(units.rad),
    n_pixels)


sigma_x = 0.5 # in units of arcsec
image = func_1D(
    x=x_rad,
    sigma_x=sigma_x * units.arcsec.to(units.rad),
)
image_FFT = np.fft.fftshift(np.fft.fft(np.fft.fftshift(image)))
u_grid = np.fft.fftshift(np.fft.fftfreq(n_pixels, d=pixel_scale * units.arcsec.to(units.rad)))

image_dFT = np.zeros(shape=n_pixels, dtype="complex")
for i in range(u_grid.shape[0]):
    for j in range(n_pixels):
        image_dFT[i] += image[j] * np.exp(
            -2.0
            * np.pi
            * 1j
            * (u_grid[i] * x_rad[j])
        )

value = 0.23

figure, axes = plt.subplots(nrows=1,ncols=3,figsize=(14,6))
axes[0].plot(x_rad * 10**6.0, image, marker="o")
for x_i in x_rad:
    axes[0].axvline(x_i * 10**6.0, linestyle="--", color="black")
axes[0].set_xlabel(r"x ($\times 10^6$; rad)")
axes[0].set_title("x-plane")

for u_grid_i in u_grid:
    axes[1].axvline(u_grid_i / 10**6.0, linestyle="--", color="black")
axes[1].plot(u_grid / 10**6.0, image_FFT.real, color="b")
axes[1].plot(u_grid / 10**6.0, image_dFT.real, color="r", linestyle="None", marker="o")
axes[1].set_title("u-plane (real)")
axes[1].set_xlabel(r"u ($\times 10^{-6}$; rad$^{-1}$)")
axes[1].plot(u_grid / 10**6.0, image_FFT.real - image_dFT.real, color="black", label="difference")

axes[2].plot(u_grid / 10**6.0, image_FFT.imag, color="b")
axes[2].plot(u_grid / 10**6.0, image_dFT.imag, color="r", linestyle="None", marker="o")
axes[2].set_title("u-plane (imag)")
axes[2].set_xlabel(r"u ($\times 10^{-6}$; rad$^{-1}$)")
#axes[2].plot(u_grid / 10**6.0, image_FFT.imag - image_dFT.imag, color="black", label="difference")
axes[1].legend()
plt.show()

Answer 1

我做了一个最小的例子（我希望）。对于 FFT 和朴素傅立叶积分（针对相同的频率值计算），我得到了基本相同的数字。

import numpy as np
import matplotlib.pyplot as p 
%matplotlib inline

def signal(x, sigma_x):
    return np.exp(-(x**2.0 / (2.0 * sigma_x**2)))

t=np.linspace(-10,10,1000)
sigma=.3
sig=np.exp(-(t**2.0 / (2.0 * sigma **2)))

p.subplot(311)
p.plot(t,sig);

ft=np.fft.fftshift(np.fft.fft(sig))
freq=np.fft.fftshift(np.fft.fftfreq(1000,0.02))
p.subplot(312)
p.plot(freq,np.abs(ft))
print(np.abs(ft)[500:505])
# naive fourier integral 
fi=[]
for f in freq: 
  i=np.sum( sig* np.exp(- 1j* 2 *np.pi*f*t ))
  fi.append(np.abs(i))

p.subplot(313)
p.plot(freq,fi)

print(np.abs(fi)[500:505])

Answer 2

我更新了@roadrunner66 的示例，改为显示 FT 的实部和虚部，而不是幅度，因为我要使用它的应用程序涉及处理 FT 的实部和虚部（通常称为作为干涉测量中的可见性）。

以下是稍微更新的示例。

import numpy as np
import matplotlib.pyplot as plt

t=np.linspace(-10,10,1000)
sigma=.3
sig=np.exp(-(t**2.0 / (2.0 * sigma **2)))

ft=np.fft.fftshift(np.fft.fft(sig))
freq=np.fft.fftshift(np.fft.fftfreq(len(t),abs(t[0] - t[1])))

# naive fourier integral
fi_real=[]
fi_imag=[]
for f in freq:
  i=np.sum( sig* np.exp(- 1j* 2 *np.pi*f*t ))
  fi_real.append(i.real)
  fi_imag.append(i.imag)

figure, axes = plt.subplots(nrows=1,ncols=2)
axes[0].plot(freq,ft.real, color="b", label="np.fft.fft")
axes[0].plot(freq,fi_real, color="r", label="exact")
axes[0].set_xlim(-5.0, 5.0)
axes[0].set_title("real")
axes[0].legend()
axes[1].plot(freq,ft.imag, color="b", label="np.fft.fft")
axes[1].plot(freq,fi_imag, color="r", label="exact")
axes[1].set_xlim(-5.0, 5.0)
axes[1].set_title("imag")
axes[1].legend()
plt.show()

查看输出图，我认为很明显 np.fft 不适合使用 FFT 的实部和虚部。