☰
光学衍射MATLAB仿真:菲涅耳、夫琅禾费与角谱法实现及避坑指南
2026/10/15 5:06:13 网站建设 项目流程

简介:这份资源面向光学、物理及光电专业的学生与教师,以及需要借助编程验证衍射理论的工程人员,提供用MATLAB模拟光学衍射实验的完整文档。压缩包内共1个doc文件,约247KB,以图文结合的方式系统整理了两类模拟思路:一是衍射积分方法,涵盖单缝衍射、杨氏双缝干涉与黑白光栅衍射的光强分布计算;二是傅立叶变换方法,涉及简单孔径、黑白光栅及正弦光栅的夫琅和费衍射模拟。文档中配有可直接运行的MATLAB代码、参数说明与衍射图样,读者能据此复现实验现象、观察干涉条纹受衍射因子调制及缺级规律,并理解屏函数功率谱与衍射场强度分布之间的对应关系。目前已有791人学习,适合作为课程实验辅助材料或自学参考。

1. 光学衍射matlab:从“公式能推”到“图能跑出来”的那道坎

很多人第一次用光学衍射matlab,卡住的地方不是公式,而是“我明明按书上的公式写了,怎么出来的图跟预期完全对不上”。这其实不是个例。光学衍射本身涉及傅里叶光学、角谱传播、菲涅耳近似、夫琅禾费近似这几套体系,每套都有自己的适用边界,而matlab里做数值计算时,采样、补零、坐标缩放、fftshift的位置,任何一个环节出错,结果都会翻车。这个方向适合两类人:一类是光学、物理、光电专业的同学,需要把课堂上的衍射积分变成可运行的仿真;另一类是做计算成像、全息显示、光束整形的工程师,需要用matlab快速验证一个传播模型是否成立。核心诉求就一个:给定波长、孔径、传播距离,能稳定地算出正确的衍射场分布,而不是每次都在调参数里碰运气。

2. 先分清三种传播模型:菲涅耳、夫琅禾费、角谱各管哪一段

2.1 从衍射积分到数值传播:为什么不能一个公式打天下

光学衍射的出发点是基尔霍夫衍射公式,它是一个面积分,直接数值算二维积分在matlab里非常慢,而且采样一密就爆内存。所以实际工程里几乎不会去硬算原始积分,而是根据传播距离和孔径尺寸的关系,切换到不同的近似模型。

判断用哪个模型,核心看一个无量纲数——菲涅耳数 (N_F = a^2/(\lambda z)),其中 (a) 是孔径特征尺寸,(\lambda) 是波长,(z) 是传播距离。当 (N_F \gg 1),也就是传播距离相对较近,用菲涅耳衍射;当 (N_F \ll 1),传播距离很远,进入夫琅禾费区,衍射场就是孔径的傅里叶变换;当传播距离极短或者你需要严格处理大角度传播,菲涅耳近似本身都不够用,这时候上角谱法。

我一般会先算一下菲涅耳数,再决定用哪套代码。很多人一上来就写菲涅耳积分,结果传播距离设成几米、孔径又是毫米级,其实早就进了夫琅禾费区,算出来的相位因子完全是多余的,图看着“差不多”但相位全错。

2.2 菲涅耳衍射的matlab实现:单步还是分步

菲涅耳衍射的数值实现有两种常见路子:单步傅里叶变换法和卷积法。单步法适合传播距离不太长、采样要求不极端的情况,代码短,跑得快。卷积法(也叫直接法)适合近距离传播,因为它对采样间隔的限制更宽松。

下面是一个单步菲涅耳衍射的matlab核心代码,假设孔径是一个矩形孔,波长532nm,传播距离0.5m:

% 菲涅耳衍射 - 单步傅里叶变换法 lambda = 532e-9; % 波长 532nm k = 2*pi/lambda; % 波数 z = 0.5; % 传播距离 0.5m L = 5e-3; % 孔径所在平面尺寸 5mm N = 1024; % 采样点数 dx = L/N; % 孔径平面采样间隔 x = (-N/2:N/2-1)*dx; [X, Y] = meshgrid(x, x); % 矩形孔径,边长1mm a = 1e-3; U0 = double(abs(X) <= a/2 & abs(Y) <= a/2); % 菲涅耳单步变换 fx = (-N/2:N/2-1)/(N*dx); % 频率坐标 [FX, FY] = meshgrid(fx, fx); H = exp(1j*k*z)/(1j*lambda*z) * exp(1j*pi*lambda*z*(FX.^2 + FY.^2)); U1 = fftshift(fft2(ifftshift(U0))) .* H; U1 = ifftshift(ifft2(fftshift(U1))); % 观察平面坐标缩放 dx1 = lambda*z/(N*dx); x1 = (-N/2:N/2-1)*dx1; I1 = abs(U1).^2; imagesc(x1*1e3, x1*1e3, I1); xlabel('x (mm)'); ylabel('y (mm)'); title('菲涅耳衍射光强分布'); axis square; colormap('hot');

这段代码里最关键的是传递函数H的构造和fftshift/ifftshift的配对。H里的二次相位因子是菲涅耳近似的核心,少了它就成了夫琅禾费衍射。dx1是观察面的采样间隔,它和孔径面采样间隔dx之间是倒数关系,这个缩放关系如果搞反,光斑尺寸会差好几个量级。

参数方面,N取1024是速度和精度的折中,如果你要观察高频细节,可以加到2048,但内存占用会翻四倍。L的选择要保证孔径被充分采样,一般让孔径占L的1/5到1/2比较合适。z不能太小,否则菲涅耳近似本身不成立,一般要求 (z > a^2/\lambda) 的若干倍。

2.3 夫琅禾费衍射:什么时候可以直接用fft2

当传播距离满足 (z \gg a^2/\lambda) 时,菲涅耳衍射公式里的二次相位因子在孔径范围内变化很小,可以忽略,衍射场直接正比于孔径透射函数的傅里叶变换。这时候代码可以简化成:

% 夫琅禾费衍射 - 直接傅里叶变换 U0 = double(abs(X) <= a/2 & abs(Y) <= a/2); U_far = fftshift(fft2(ifftshift(U0))); I_far = abs(U_far).^2; imagesc(I_far); colormap('hot'); title('夫琅禾费衍射图样');

注意这里没有乘任何相位因子,因为夫琅禾费区我们通常只关心强度分布,相位因子在强度里被消掉了。但如果你要做干涉或者全息,相位不能丢,那就还是得回到菲涅耳或者角谱。

夫琅禾费衍射的观察面坐标和孔径面坐标之间是傅里叶变换对的关系,实际物理尺寸需要根据焦距或者传播距离来换算。很多人直接拿fft2的结果当物理坐标用,这是不对的,fft2输出的索引对应的是空间频率,要乘以 (\lambda z) 才能变成实际坐标。

2.4 角谱法:宽角度、短距离传播的可靠选择

角谱法不依赖菲涅耳近似,它把光场分解成不同方向传播的平面波,每个平面波乘以对应的传播相位因子,再叠加回去。角谱法的传递函数是:

[ H(f_x, f_y) = \exp\left(j 2\pi z \sqrt{1/\lambda^2 - f_x^2 - f_y^2}\right) ]

当 (f_x^2 + f_y^2 > 1/\lambda^2) 时,对应倏逝波,传播因子变成衰减的,数值上要处理成零或者衰减项。

% 角谱法传播 lambda = 532e-9; k = 2*pi/lambda; z = 0.01; % 近距离 10mm L = 2e-3; N = 1024; dx = L/N; x = (-N/2:N/2-1)*dx; [X, Y] = meshgrid(x, x); U0 = double(abs(X) <= 0.5e-3 & abs(Y) <= 0.5e-3); fx = (-N/2:N/2-1)/(N*dx); [FX, FY] = meshgrid(fx, fx); arg = 1/lambda^2 - FX.^2 - FY.^2; H = exp(1j*2*pi*z*sqrt(max(arg, 0))); H(arg < 0) = 0; % 倏逝波截止 U1 = ifft2(ifftshift(fftshift(fft2(ifftshift(U0))) .* H)); I1 = abs(U1).^2; imagesc(x*1e3, x*1e3, I1); colormap('hot');

角谱法的采样要求比菲涅耳严格,dx必须小于 (\lambda/2) 才能避免混叠,但实际中为了兼顾视场,往往取dx在 (\lambda/2) 到 (\lambda) 之间,这时候高频部分会有误差,但低频结构还是准的。z很小时角谱法优势明显,菲涅耳近似在z小于孔径尺寸时误差会很大。

3. 采样、补零与坐标缩放:三个让结果对不上的隐形杀手

3.1 采样间隔怎么定:从孔径尺寸反推N

采样间隔dx不是随便设的。它必须满足两个条件:第一,孔径的精细结构要被采到,一般要求dx小于最小特征尺寸的1/2;第二,观察面的采样间隔dx1 = λz/(Ndx)要能覆盖你关心的光斑范围。

这两个条件有时候是矛盾的。比如孔径很小、传播距离很远,dx1会很大,光斑在观察面上只占几个像素,看着很粗糙。这时候要么增大N,要么在孔径面补零。补零相当于在频域插值,能让观察面的采样变密,但不会增加真实信息量。

我一般会先算一下观察面的期望光斑尺寸,然后反推需要的dx1,再根据dx1 = λz/(Ndx)确定N。如果N太大跑不动,就考虑用两步法或者缩放傅里叶变换。

3.2 补零的正确姿势:位置比数量更重要

补零在matlab里做fft时很常见,但补零的位置会影响结果的相位。如果你用ifftshift把零频移到角落再补零,补完再fftshift回来,相位是对的。如果直接在矩阵中间插零,相位就乱了。

% 正确的补零方式 N0 = 512; Npad = 2048; U0_pad = zeros(Npad, Npad); start = (Npad - N0)/2 + 1; U0_pad(start:start+N0-1, start:start+N0-1) = U0;

补零之后,频率坐标也要相应更新,dx不变但N变了,频率间隔变成1/(Npad*dx),传递函数H要用新的频率坐标重新算。很多人补零之后忘了改H,结果衍射图样完全不对。

3.3 坐标缩放:观察面的物理尺寸到底是多少

这是最容易出错的地方。fft2输出的矩阵索引是整数,对应的物理坐标是 (x_m = m \cdot \Delta x_1),其中 (\Delta x_1 = \lambda z/(N \Delta x))。这个公式只对菲涅耳单步法和夫琅禾费衍射成立,角谱法的观察面坐标和孔径面坐标是一样的,不需要缩放。

如果你用的是卷积法或者两步法,坐标缩放关系又不一样。我见过有人把菲涅耳的缩放公式套到角谱法上,结果光斑尺寸差了100倍,还以为是衍射效应。

提示:每次换传播模型,先拿一个已知解析解的情况验证坐标。比如圆孔夫琅禾费衍射的第一暗环半径是 (1.22\lambda z/D),算出来和仿真对一下,对上了再往下做。

4. 避坑与排查:光学衍射仿真里最常见的五个翻车现场

4.1 光斑出现十字条纹或者棋盘格

现象:衍射图样上叠加了规则的十字或者网格状条纹,不是物理衍射该有的样子。

原因:fftshift和ifftshift用反了,或者该用的时候没用。matlab的fft2默认零频在左上角,而物理上我们习惯零频在中心。如果传递函数H是在中心化的频率坐标下算的,但fft2没有先ifftshift,H和频谱就对不上,产生高频干涉条纹。

解决:记住一个原则——fft2之前用ifftshift把零频移到左上角,fft2之后用fftshift把零频移回中心。传递函数H如果在中心化坐标下定义,要先ifftshift再乘。代码里养成成对出现的习惯,不要凭感觉。

4.2 传播距离增大后光斑反而变小

现象:按菲涅耳衍射,传播距离越大光斑应该越大,但仿真里z增大后光斑缩小了。

原因:坐标缩放公式用反了。菲涅耳单步法的观察面采样间隔是 (\Delta x_1 = \lambda z/(N\Delta x)),z增大时(\Delta x_1)增大,但如果你用imagesc画图时坐标轴还是用的孔径面坐标,视觉上光斑就“变小”了。实际上光斑物理尺寸是变大的,只是坐标轴没更新。

解决:画图时一定要用观察面的实际物理坐标,不要图省事直接用像素索引。把x1算出来,imagesc(x1, x1, I1),这样z变化时坐标轴会自动缩放。

4.3 高频细节丢失或者出现混叠

现象:孔径的锐利边缘在衍射图样里变得模糊,或者出现不该有的高频波纹。

原因:采样不足。孔径面的dx太大,孔径边缘没有被充分采样,相当于给孔径加了一个低通滤波。或者观察面的dx1太大,衍射图样的高频部分被截断。

解决:先检查dx是否小于孔径最小特征的1/2。如果孔径有锐利边缘,考虑用超高斯或者平滑窗代替硬边,减少高频分量。如果观察面采样不够,增大N或者用补零。混叠的话,检查频率坐标是否超出了奈奎斯特频率。

4.4 角谱法在远距离传播时结果发散

现象:用角谱法算远距离传播,光斑越来越大但强度不守恒,甚至出现数值溢出。

原因:角谱法的传递函数是纯相位因子,模为1,理论上能量守恒。但如果频率坐标超出了倏逝波截止频率没有置零,sqrt里出现负数,matlab会返回复数,导致结果异常。另外远距离传播时,角谱法的采样要求变得很苛刻,dx需要非常小,否则高频相位变化太快,采样跟不上。

解决:确保arg < 0的位置H置零。远距离传播优先用菲涅耳或者夫琅禾费,角谱法适合近距离。如果非要用角谱法传远距离,考虑分段传播,每段距离短一些,中间可以重采样。

4.5 相位分布看起来像噪声

现象:强度图看着正常,但相位图完全是随机的,没有物理意义。

原因:强度为零的区域,相位是没有定义的,matlab算出来的相位就是数值噪声。另外如果用了fft2而没有正确处理相位因子,相位会带上一个线性或者二次的额外项。

解决:看相位时把强度低于阈值的区域掩掉,只显示强度足够高的区域的相位。检查传递函数的相位符号是否正确,菲涅耳传递函数的二次相位因子符号很容易搞反,正负号错了相位就完全不对。

5. 从能跑到好用:一个可复用的衍射仿真函数框架

5.1 把三种模型封装成一个函数

每次写脚本都复制粘贴容易出错,我一般会封装一个函数,输入波长、孔径、传播距离、模型类型,输出观察面的光场和坐标。

function [U1, x1, y1] = diffraction_propagate(U0, dx, lambda, z, model) % U0: 输入光场 % dx: 孔径面采样间隔 % lambda: 波长 % z: 传播距离 % model: 'fresnel', 'fraunhofer', 'angular' [N, ~] = size(U0); k = 2*pi/lambda; fx = (-N/2:N/2-1)/(N*dx); [FX, FY] = meshgrid(fx, fx); switch model case 'fresnel' H = exp(1j*k*z)/(1j*lambda*z) * ... exp(1j*pi*lambda*z*(FX.^2 + FY.^2)); U1 = ifftshift(ifft2(fftshift( ... fft2(ifftshift(U0)) .* ifftshift(H)))); dx1 = lambda*z/(N*dx); x1 = (-N/2:N/2-1)*dx1; case 'fraunhofer' U1 = fftshift(fft2(ifftshift(U0))); dx1 = lambda*z/(N*dx); x1 = (-N/2:N/2-1)*dx1; case 'angular' arg = 1/lambda^2 - FX.^2 - FY.^2; H = exp(1j*2*pi*z*sqrt(max(arg, 0))); H(arg < 0) = 0; U1 = ifft2(ifftshift(fftshift(fft2(ifftshift(U0))) .* H)); x1 = (-N/2:N/2-1)*dx; end y1 = x1; end

这个框架的好处是模型切换只改一个参数,坐标缩放自动处理。注意角谱法的坐标不缩放,菲涅耳和夫琅禾费缩放,这个差异在函数里已经处理了。

5.2 用解析解验证你的代码

写完函数不要急着跑复杂孔径,先用圆孔验证。圆孔夫琅禾费衍射的强度分布有解析解,第一暗环半径是 (r_1 = 1.22\lambda z/D)。仿真跑出来,找到第一暗环的位置,和解析解对比,误差在5%以内说明采样和坐标缩放都对了。

% 圆孔验证 D = 1e-3; U0 = double(sqrt(X.^2 + Y.^2) <= D/2); [U1, x1, ~] = diffraction_propagate(U0, dx, lambda, z, 'fraunhofer'); I1 = abs(U1).^2; I1 = I1 / max(I1(:)); % 找第一暗环 center = N/2 + 1; profile = I1(center, center:end); [~, idx] = min(profile(1:round(N/4))); r_sim = x1(center + idx - 1); r_theory = 1.22*lambda*z/D; fprintf('仿真第一暗环半径: %.4f mm\n', r_sim*1e3); fprintf('理论第一暗环半径: %.4f mm\n', r_theory*1e3);

如果误差大,先检查dx是否足够小,再检查z是否真的满足夫琅禾费条件。有时候z不够大,菲涅耳效应还没完全消失,第一暗环位置会偏移,这是物理上的偏差,不是代码错了。

5.3 参数扫描和批量仿真

实际工作中经常需要扫传播距离或者扫孔径尺寸,看衍射图样怎么变化。这时候把上面的函数放进循环,注意每次循环重新算H,因为H和z有关。

z_list = [0.1, 0.5, 1, 2]; figure; for i = 1:length(z_list) [U1, x1, ~] = diffraction_propagate(U0, dx, lambda, z_list(i), 'fresnel'); subplot(2, 2, i); imagesc(x1*1e3, x1*1e3, abs(U1).^2); colormap('hot'); axis square; title(sprintf('z = %.1f m', z_list(i))); end

批量仿真时注意内存,N=2048时一个复数矩阵就是2048×2048×16字节,约64MB,循环里及时clear不用的变量。如果z_list很长,考虑把结果存成视频或者gif,不要全部留在内存里。

5.4 一个容易忽略的细节:光场归一化

衍射传播不改变总能量,但数值上fft2会引入缩放因子。如果你关心绝对强度,需要在传播后做归一化,让总能量和输入一致。如果只关心相对分布,归一化可以不做,但画图时最好除以最大值,避免不同z之间没法对比。

% 能量归一化 U1 = U1 * sqrt(sum(abs(U0(:)).^2) / sum(abs(U1(:)).^2));

这个操作在比较不同传播距离的衍射效率时很有用。我见过有人忘了归一化,结果远距离的衍射图样看着比近距离暗很多,以为是能量损失,其实是fft2的缩放因子在作怪。

光学衍射matlab仿真这件事,说到底就是三件事:选对模型、配对坐标、验证解析解。模型选错了后面全白搭,坐标配错了图看着像那么回事但数值全错,不验证解析解你永远不知道自己对不对。我自己的习惯是每写一个新模型,先拿圆孔或者单缝跑一遍,和教科书上的图对一下,对上了再往复杂孔径上套。这个习惯帮我省了很多后悔药。希望帮到你。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询