5

我写了一个简短的 matlab 脚本文件，假设运行菲涅耳的传播（衍射），这样给定某个输入字段 U0，它会告诉您该字段在距离 z0 后的样子。我将结果与教科书的结果进行了比较，看来我的程序运行良好。问题是如果我尝试采取两个传播步骤而不是一个。即，我没有对程序进行一次迭代来传播距离 z0，而是对程序进行两次迭代来传播距离 z0/2。然后我完全胡说八道，我无法弄清楚问题出在哪里。任何建议都将被非常感激地接受。这是代码：

function U = fresnel_advance (U0, dx, dy, z, lambda)
% The function receives a field U0 at wavelength lambda
% and returns the field U after distance z, using the Fresnel
% approximation. dx, dy, are spatial resolution.

k=2*pi/lambda;
[ny, nx] = size(U0); 

Lx = dx * nx;
Ly = dy * ny;

dfx = 1./Lx;
dfy = 1./Ly;

u = ones(nx,1)*((1:nx)-nx/2)*dfx;    
v = ((1:ny)-ny/2)'*ones(1,ny)*dfy;   

O = fftshift(fft2(U0));

H = exp(1i*k*z).*exp(-1i*pi*lambda*z*(u.^2+v.^2));  

U = ifft2(O.*H);
4

4 回答 4

5

呼叫后fft2，您呼叫也fftshift使直流频率在中间。

但是当您调用 时ifft2，该函数假定您的直流频率仍为 (1,1)。因此，在进行逆 FFT 之前，您必须回到这种格式。

因此将最后一行更改为U = ifft2(fftshift(O.*H))可能会解决问题。

编辑

我刚刚看到 Matlab 建议在两次 insteafifftshift之后使用（找不到引入它的版本）。根据文档，在奇数大小的情况下，调用序列和不等价。fftshiftifftshiftifftshift(fftshift(X))ifftshift(fftshift(X))

所以我认为最好这样做：U = ifft2(ifftshift(O.*H))在代码的最后一步。

于 2014-01-07T13:40:45.503 回答
2

如果您可以将一些示例输入发布到您的程序中以演示该问题，那将会很有帮助。

我怀疑这个问题可能与您调用 FFTSHIFT 的次数不够多有关。通常，人们认为光场矩阵的中心位于“原点”上，而 FFT2 则认为是“左下角”。因此，您应该在 FFT2 之前和之后进行 FFTSHIFT。

你也需要为 IFFT2 做同样的事情。

编辑 向 FFTSHIFT 添加两个调用的理由：比较这两个：

N = 512; [x,y] = meshgrid(-1:1/N:(N-1)/N);
mask = (x.*x + y.*y) < 0.001;
figure(1)
imagesc(angle(fftshift(fft2(fftshift(mask)))))
figure(2)
imagesc(angle(fftshift(fft2(mask)))
于 2014-01-07T13:13:44.760 回答
2

事实上，问题在于您运行 fft 的方式。Khedar Kare, Wiley 2015在傅里叶光学和计算成像中对此进行了很好的解释：

因此，在大多数编程平台中，为了使结果从物理角度有意义（例如在描述衍射等现象时），2D FFT 的适当序列由以下公式给出：fftshift(fft2(ifftshift(...)))。

在您的代码中，您应该：O = fftshift(fft2(ifftshift(U0)));

如果您对用 Python 开发的软件感兴趣，这里有一个用于光学（包括衍射）的快速开发的 Python 工具箱：PyOptica。在 PyOptica 中，可以使用以下方法传播波前：

import astropy.units as u
import numpy as np

import pyoptica as po


wavelength = 500 * u.nm 
pixel_scale = 22 * u.um
npix = 1024

w = 6 * u.mm
h = 3 * u.mm
axis_unit = u.mm

wf = po.Wavefront(wavelength, pixel_scale, npix)
ap = po.RectangularAperture(w, h)
wf = wf * ap
fig_1 = po.plotting.plot_wavefront(wf, 'intensity', axis_unit=axis_unit)

初始孔径

第二步是传播：

f = 50 * u.cm
fs_f = po.FreeSpace(f)
wf_forward = wf * fs_f
fig_2 = po.plotting.plot_wavefront(wf_forward, 'intensity',  axis_unit=axis_unit)

传播后的结果

记住菲涅耳传播的采样条件很重要： 在此处输入图像描述( z <= N (dx)^2 / lambda) 其中：

  • N——像素数（一个方向）；
  • dx——像素大小；
  • λ - 波长。

该条件基于计算傅里叶光学： David Voelz 的 MATLAB 教程，SPIE 2011

您应该在代码中实现该条件。在 PyOptica 中，它总是在传播之前进行检查；如果请求的距离违反条件，则传播距离被分解为子步。

于 2020-03-22T13:03:18.467 回答
0

我认为相位项有错误。"H = exp(1i k z) .exp(-1i pi lambda z*(u.^2+v.^2)); "
应该是 "H = exp(1i k z) .exp(-1i pi/ λ/z*(u.^2+v.^2)); "

于 2018-04-25T17:11:59.163 回答