在快速傅立叶变换的频率下计算dFT

Question

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

我是不是遗漏了什么？根据定义，它们不应该是相同的吗？(当我查看dFT/FFT的图像部分时，差别更大)。

我附上了我使用的脚本，它生成了这个比较dFT和快速傅立叶变换的真实和图像部分的图。

​

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

我通过@ real runner66更新了示例，以显示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是不合适的。

​

​