【问题标题】:2d Fourier Transforms: FFT vs Fourier Optics二维傅里叶变换:FFT 与傅里叶光学
【发布时间】:2021-07-12 21:26:15
【问题描述】:

我正在尝试使用编程来增加我对傅里叶光学的理解。我知道傅里叶变换的傅里叶变换在物理上和数学上是反转的 -> F{F{f(x)} = f(-x)。我有两个问题1)第二个变换不会返回任何像原始函数一样的东西,除了简单的高斯情况(这使得它更加混乱),并且 2)似乎有一些缩放需要我“放大”并将转换后的图像扭曲到一个没有帮助的程度的因素(如下图所示)。 **根据@Cris Luengo 的建议编辑

#%% Playing with 2d Fouier Transform

import numpy as np
from scipy import fftpack
import matplotlib.pyplot as plt
import LightPipes

wavelength = 792*nm
size = 15*mm
N = 600
w0=3*mm

# Fields
sq = np.zeros([100,100])
sq[25:75, 25:75] = 1
F=Begin(size,wavelength,N)
I0 = Intensity(0,GaussBeam(F, w0, LG=True, n=0, m=0))
I1 = Intensity(0,GaussBeam(F, w0, LG=False, n=0, m=1))+Intensity(0,GaussBeam(F, w0, LG=False, n=1, m=0))

# Plot transforms
f = sq
F = np.fft.fftshift(fftpack.fft2(f))
F_F = fftpack.fft2((F))

plt.subplot(331), plt.imshow(f)
plt.title(r'f'), plt.xticks([]), plt.yticks([])
plt.subplot(332), plt.imshow(np.abs(F))
plt.title(r'F\{f\}'), plt.xticks([]), plt.yticks([])
plt.subplot(333), plt.imshow(np.abs(F_F))
plt.title('F\{F\{f\}\}'), plt.xticks([]), plt.yticks([])

# plt.subplot(331), plt.imshow(f)
# plt.title(r'f'), plt.xticks([]), plt.yticks([])
# plt.subplot(332), plt.imshow(np.abs(F))
# plt.title(r'F\{f\}'), plt.xticks([]), plt.yticks([])
# plt.subplot(333), plt.imshow(np.abs(F_F))
# plt.title('F\{F\{f\}\}'), plt.xticks([]), plt.yticks([])

f = I0
F = np.fft.fftshift(fftpack.fft2(f))
F_F = fftpack.fft2((F))

plt.subplot(334), plt.imshow(f)
plt.title(r'f'), plt.xticks([]), plt.yticks([])
plt.subplot(335), plt.imshow(np.abs(F))
plt.title(r'F\{f\}'), plt.xticks([]), plt.yticks([])
plt.subplot(336), plt.imshow(np.abs(F_F))
plt.title('F\{F\{f\}\}'), plt.xticks([]), plt.yticks([])

f = I1
F = fftpack.fft2(f)
F_F = fftpack.fft2(F)

plt.subplot(337), plt.imshow(f)
plt.title(r'f'), plt.xticks([]), plt.yticks([])
plt.subplot(338), plt.imshow(np.abs(F))
plt.title(r'F\{f\}'), plt.xticks([]), plt.yticks([])
plt.subplot(339), plt.imshow(np.abs(F_F))
plt.title('F\{F\{f\}\}'), plt.xticks([]), plt.yticks([])

plt.tight_layout()
plt.show()

【问题讨论】:

  • 如果您删除对fftshiftabs 的调用,事情应该可以正常工作。 I0I1 是什么?您的代码中是否需要这么多冗余来说明您的问题?
  • @Cris Luengo 抱歉,我在为此制作的演示中发现了一个错误,并且不小心删除了显示 I0I1 是什么的图片。我相信fftshift 是正确重新排列象限所必需的。关于标准的东西是将高频放在边缘?并且np.abs 是必要的,因为fft2 的结果很复杂。
  • 我试图礼貌地暗示我一开始就尝试过。由于这个问题,abs 在那里。在 stackoverflow TypeError: Image data of dtype complex128 cannot be converted to float 上找到了对它的更正。这种转变来自 youtube 对 FFT 的描述。在发现当使用小(10x10)正方形时,强度在角落是错误的。它们必须被转移。
  • @CrisLuengo 这是促使我尝试添加np.abs stackoverflow.com/a/38333442/14305494的答案
  • fftshift 是将原点从左上角(DFT/FFT 期望的位置)移动到我们喜欢看到它的中心。仅当您想要显示 FFT 的结果时才使用它。 abs 丢弃 DFT 的相位,破坏您的数据。它会导致所有正弦分量在原点对齐,从而导致每个结果中的特征单峰。 fft(fft(f)) 将产生您期望的结果。它不会完全是f(-x),因为我们处理的是离散傅里叶变换,而不是普通的 FT。

标签: python scipy fft


【解决方案1】:

和 Cris 聊天后,好像没有缩放因子,这种类型的 DFT 似乎就是这样工作的。所以我找到的解决方案是将像素增加到可以放大并获得足够清晰图像的程度。这不是一个很好的解决方案,但与LightPipes 搭配使用,现在可以了解光模式的变换会是什么样子,并说明在镜头系统的图像平面上它们会像在前焦场。

#%% Playing with 2d Fourier Transform

import numpy as np
import matplotlib.pyplot as plt
import LightPipes
from scipy.fftpack import fft2 as fft
from numpy.fft import fftshift, ifftshift

wavelength = 792*nm
size = 100*mm
N = 1000
w0=3*mm

# Fields
sq = np.zeros([100,100])
sq[25:75, 25:75] = 1
F=Begin(size,wavelength,N)
I0 = Intensity(0,GaussBeam(F, w0, LG=True, n=0, m=0))
I1 = Intensity(0,GaussBeam(F, w0, LG=False, n=0, m=1))+Intensity(0,GaussBeam(F, w0, LG=False, n=1, m=0))

# Plot transforms
f = sq
F = fftshift(fft(ifftshift(f)))
F_F = fftshift(fft(ifftshift(F)))

plt.subplot(331), plt.imshow(f,cmap = cmap)
plt.title(r'f'), plt.xticks([]), plt.yticks([])
plt.subplot(332), plt.imshow(np.abs(F),cmap = cmap)
plt.title(r'F\{f\}'), plt.xticks([]), plt.yticks([])
plt.subplot(333), plt.imshow(np.abs(F_F),cmap = cmap)
plt.title('F\{F\{f\}\}'), plt.xticks([]), plt.yticks([])

f = I0
F = fftshift(fft(ifftshift(f)))
F_F = fftshift(fft(ifftshift(F)))

plt.subplot(334), plt.imshow(f[450:550,450:550],cmap = cmap)
plt.title(r'f'), plt.xticks([]), plt.yticks([])
plt.subplot(335), plt.imshow(np.abs(F)[450:550,450:550],cmap = cmap)
plt.title(r'F\{f\}'), plt.xticks([]), plt.yticks([])
plt.subplot(336), plt.imshow(np.abs(F_F)[450:550,450:550],cmap = cmap)
plt.title('F\{F\{f\}\}'), plt.xticks([]), plt.yticks([])

f = I1
F = fftshift(fft(ifftshift(f)))
F_F = fftshift(fft(ifftshift(F)))

plt.subplot(337), plt.imshow(f[450:550,450:550],cmap = cmap)
plt.title(r'f'), plt.xticks([]), plt.yticks([])
plt.subplot(338), plt.imshow(np.abs(F)[450:550,450:550],cmap = cmap)
plt.title(r'F\{f\}'), plt.xticks([]), plt.yticks([])
plt.subplot(339), plt.imshow(np.abs(F_F)[450:550,450:550],cmap = cmap)
plt.title('F\{F\{f\}\}'), plt.xticks([]), plt.yticks([])

plt.tight_layout()
plt.show()

【讨论】:

    猜你喜欢
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2012-05-18
    • 1970-01-01
    • 2019-11-10
    相关资源
    最近更新 更多