基于FFT的夫琅禾费衍射与光栅MATLAB仿真统一框架
2026/9/16 19:01:07 网站建设 项目流程

简介:一套基于MATLAB的光学仿真资源,围绕夫郎费禾衍射、双缝干涉、平面光栅衍射和单缝衍射四个经典物理光学现象展开,适合本科、硕士阶段教研学习,也可作为课程设计与毕业设计的参考实现。压缩包共12个文件,包含4个MATLAB脚本(.m)、7张运行结果图(.png)以及1个说明文档(.txt),整体大小仅484KB,轻量易用。脚本封装了从参数设置、数值计算到图像绘制的完整流程,运行后可直接得到衍射与干涉图样;通过调整缝宽、波长等参数,能快速观察不同条件下的条纹变化与光强分布,配合截图和说明文档,便于对照理论公式加深理解。已有258人学习/下载,比较适合正在学习光学仿真、希望快速上手MATLAB数值模拟的读者。

1. 夫琅禾费衍射、双缝、平面光栅和单缝,为什么都用同一个 MATLAB 仿真套路

夫琅禾费衍射、双缝干涉、平面光栅和单缝这四类经典光学现象,教材里都能给出解析公式:sinc 平方、双光束干涉调制、光栅方程。但真要用 MATLAB 做光学仿真,我不建议只照着公式去画曲线。原因很直接:公式套图的代码只能验证某个特定参数组合,换一个孔径形状、换一种入射条件就要重推数学;而以傅里叶变换为核心的仿真流程,一次建模就能覆盖单缝、双缝、平面光栅,甚至任意二维掩膜。这里最该建立的一个反直觉认知是:远场夫琅禾费衍射图样就是孔径复振幅的傅里叶变换。所以这篇博客把重点放在三件事上:如何正确构造孔径掩膜,如何把 FFT 输出的频率轴换算成屏幕物理坐标,以及如何用同一套参数快速排掉仿真发散、位置错位这些常见问题。适合光学实验课、MATLAB 课程设计和想用衍射来理解 FFT 的工程师。

2. 夫琅禾费衍射模型与仿真网格参数:从远场公式到 FFT 坐标映射

2.1 夫琅禾费的物理条件与傅里叶变换表示

夫琅禾费衍射成立的条件是观察距离足够远,或者用一个会聚透镜把无穷远的衍射场拉到后焦面。习惯上写成

[ E_{far}(u,v)=C\iint E_0(x,y)\exp[-j2\pi(ux+vy)],dxdy ]

其中 (u=x_s/(\lambda f)),(v=y_s/(\lambda f)),(f) 是等效焦距,(x_s,y_s) 是屏幕上物理坐标。这个式子说明,夫琅禾费衍射场的复振幅分布就是孔径函数的傅里叶变换,观测屏上的光强则是这个变换结果的模平方。MATLAB 仿真的核心不是去直接离散这个积分,而是把孔径写成数组,做一次fft,再把频率轴按 (x_s=\lambda f \cdot f_x) 映射回屏幕坐标。

这里要注意的是:如果仿真单缝、双缝和平面光栅,它们只在横向一维有结构,用一维fft就足够;但如果想仿真圆孔、方孔、光栅二维点阵,则需要把fft换成fft2,并生成二维meshgrid坐标。整体思路完全不变,工作量只多了一个维度。很多把仿真做“乱”的初学者,就是在这个点上没有区分一维与二维,导致明明是条状结构,却因为用了二维 FFT 而出现两个方向的衍射条纹,反而把自己绕晕。

2.2 离散网格上如何设定波长、焦距和屏幕坐标

离散仿真要给出四组参数:波长、液晶等效焦距、孔径面计算区域边长 (D)、采样点数 (N)。波长决定条纹尺度,焦距决定坐标缩放,而 (D) 和 (N) 共同决定屏幕能看到的范围与分辨率。

FFT 生成的频率轴满足如下对应关系。孔径面空间步长是 (dx=D/N),FFT 频率分辨率是 (1/D),最大可表示频率是 (1/(2dx))。映射到屏幕物理坐标后,屏幕像素间距为 (dx_{obs}=\lambda f/D),屏幕横坐标总范围近似为 (\pm\lambda f/(2dx))。这个关系是整个仿真最容易被忽略的地方:屏幕坐标不是随便选的,它由波长、焦距和孔径面网格共同决定。

以一组常用参数为例:

参数取值说明
波长 (\lambda)632.8e-9 mHe-Ne 激光
焦距 (f)1.0 m等效远场距离
计算区域 (D)10e-3 m孔径面宽度
采样点数 (N)4096取 2 的幂
孔径面步长 (dx)2.44e-6 m(D/N)
屏幕像素间距 (dx_{obs})6.33e-5 m(\lambda f/D)
屏幕半宽约 0.13 m(\lambda f/(2dx))

另一个不容易想到的细节是零填充。缝的实际宽度可能只有 50 μm,而计算区域 (D) 取 10 mm,这意味着孔径数组里大部分元素是 0。零填充等效于在频域做插值,让衍射图案更平滑,但并不会提高物理上的分辨能力。真正决定峰值宽度的还是孔径物理尺寸和参与干涉的缝数量。理解这一点,才能解释为什么同一个代码改大 (N) 后条纹只是变光滑,而不是变多。

2.3 用 MATLAB 复现远场衍射的最小骨架代码

lambda = 632.8e-9; % 波长 f = 1.0; % 等效焦距 D = 10e-3; % 孔径面计算区域 N = 4096; % 采样点 dx = D/N; x = (-N/2:N/2-1)*dx; % 孔径面坐标 fx = (-N/2:N/2-1)/D; % 空间频率轴 x_obs = lambda * f * fx; % 屏幕物理坐标 % 单缝孔径 a = 50e-6; E0 = double(abs(x) <= a/2); % 夫琅禾费仿真 E_far = fftshift(fft(ifftshift(E0))); I = abs(E_far).^2; % 查看结果 plot(x_obs, I/max(I)); xlabel('x_screen / m'); ylabel('归一化强度');

这里的ifftshiftfftshift必须成对使用。孔径数组是按照 (x=0) 放在数组中央的方式构造的,直接调用fft会把原点当数组第一个元素,导致结果左右颠倒;ifftshift先把原点移到 FFT 所需的初始位置,fft计算后再用fftshift移回来。频率轴fx=(-N/2:N/2-1)/D也和这个索引顺序一致。如果不先确认这一行,后面所有衍射峰位置都会对不上。

3. 单缝衍射与双缝干涉:一个 fft 脚本把两种经典图样做出来

3.1 单缝:矩形孔径的傅里叶变换与 sinc 强度

单缝是最适合验证整套参数配置的案例。缝宽 (a) 确定后,解析解是 sinc 函数平方:

[ I=I_0\left(\frac{\sin\beta}{\beta}\right)^2,\quad \beta=\frac{\pi a x_s}{\lambda f} ]

用 MATLAB 做数值仿真时,只需要把孔径掩膜设成一个符合abs(x)<=a/2的数组。

% 使用 2.3 节的 x, fx, x_obs, lambda, f, D, N a = 50e-6; E0 = double(abs(x) <= a/2); E_far = fftshift(fft(ifftshift(E0))); I_single = abs(E_far).^2; plot(x_obs, I_single/max(I_single)); xlabel('x_screen / m'); ylabel('归一化强度');

运行后会看到中央主极大两侧等间距分布着次极大,次极大强度从 4.7% 开始逐级衰减,零点的理论位置在 (x_m=m\lambda f/a)。以当前参数计算,第一个零点约在 12.66 mm 处,屏幕半宽约 129 mm,能完整看到几个旁瓣。如果画出来后发现中央主峰旁边全是高频毛刺,通常说明N太小或者孔径周围存在非零边界,也就是数组在计算区域边缘被截断,FFT 把不连续处当成新孔径重新衍射。处理办法是加宽D,让孔径函数在数组边界处衰减到 0。

3.2 双缝:两个矩形窗叠加后的干涉调制

双缝本质上是一个孔径函数,由两个矩形脉冲相加得到,两缝中心距为 (d),缝宽仍为 (a)。傅里叶变换后,复振幅正比于 (\cos(\pi d x_s/(\lambda f))) 与单缝 sinc 函数的乘积,所以光强是干涉条纹和单缝衍射包络相乘的结果。

d = 500e-6; % 缝间距,中心到中心 a = 50e-6; % 缝宽 E0 = (abs(x - d/2) <= a/2) | (abs(x + d/2) <= a/2); E0 = double(E0); E_far = fftshift(fft(ifftshift(E0))); I_double = abs(E_far).^2; plot(x_obs, I_double/max(I_double)); xlabel('x_screen / m'); ylabel('归一化强度');

运行结果中,细密的干涉条纹周期是 (\Delta x=\lambda f/d),当前参数约为 1.27 mm;而整体亮暗由单缝的 sinc 包络控制,包络零点仍在 (\lambda f/a\approx12.66) mm 附近。这个现象恰好解释了为什么双缝干涉图样中央亮纹两侧的亮度不是均匀衰减的。很多初学者把双缝直接当成两个无限细的 δ 函数,出来的条纹确实等亮,但只适用于缝宽远小于波长的极限;实际实验中,缝宽对干涉级次有明显调制作用。

3.3 参数表:缝宽、缝距、屏幕范围怎么选

想看到的图像推荐参数判断标准
单缝旁瓣清晰(a=100) μm,(N=4096)主瓣半宽约 6.33 mm
双缝条纹明显(a=50) μm,(d=500) μm条纹间距 1.27 mm,包络零点 12.66 mm
双缝的包络也可见(a=30) μm,(d=250) μm包络零点 21 mm,包含多个干涉周期
光栅多级衍射(d=50) μm,(b=25) μm级次间距 12.66 mm,屏幕半宽 129 mm

屏幕范围由 (D) 和 (N) 共同决定。如果发现双缝条纹在屏幕边缘仍很密,不要直接加大倍率去“放大”,而是应该调整 (\lambda f/d) 的比例。比如d太小,屏幕上看不清;D取得太小,又会把观察窗压缩导致只能看到中央几条。先画一次完整图,再逐项调节Dda,比盲目改N更有效率。

3.4 仿真发散或条纹错位的常规检查清单

仿真不会像迭代求解那样发散,但会出现“看起来发散”的高频噪声或完全错误的条纹结构。遇到这类问题,我一般按顺序检查四件事:

第一,ifftshiftfftshift是否成对。漏掉任何一个,图形可能左右翻转或中心偏移。第二,频率轴是否用的是(-N/2:N/2-1)/D,而不是(-N/2:N/2-1)/dx。后者会把所有峰的物理位置放大或缩小,出现“坐标对不上”。第三,孔径是否超界。比如abs(x) <= a/2本身没问题,但如果在循环里累加多个缝,要确保所有非零区域都在[-D/2, D/2)内,否则 FFT 会把数组右侧当成左侧继续周期延拓。第四,观察强度时是否用了abs(E_far).^2。如果只画real(E_far),会看到大量正负交替的旁瓣,那是相位而不是光强。

4. 平面光栅衍射仿真:把多缝叠加成衍射级,并验证光栅方程

4.1 平面透射光栅的孔径构造

平面透射光栅可以看成多个等间距单缝的集合,缝数为 (M),缝间距为 (d),每一条缝宽为 (b)。构造孔径掩膜时,只要在x坐标上循环生成多个矩形脉冲即可。

M = 9; % 缝数 d = 50e-6; % 光栅常数 b = 25e-6; % 缝宽 E0 = zeros(1, N); centers = (-(M-1)/2 : (M-1)/2) * d; for c = centers E0(abs(x - c) <= b/2) = 1; end E_far = fftshift(fft(ifftshift(E0))); I_grating = abs(E_far).^2; plot(x_obs, I_grating/max(I_grating)); xlabel('x_screen / m'); ylabel('归一化强度');

这里选择奇数的M=9,是为了让光栅中心有一条缝正好位于 (x=0),这样中央衍射级的相位中心明确。centers从 (-4d) 到 (4d),总宽度约为 400 μm,远小于计算区域 (D=10) mm,因此边界截断效应很小。b=25e-6是缝宽,和光栅常数 (d=50e-6) 的比值为 1:2,这个比例会在后面引出缺级现象。

4.2 用光栅方程核对各级衍射峰位置

平面光栅在正入射情况下的主极大条件为:

[ d\sin\theta_m=m\lambda ]

在傍轴条件下,屏幕上峰位近似为:

[ x_m\approx\frac{m\lambda f}{d} ]

当前参数下,(m=1) 对应 12.66 mm,(m=2) 对应 25.31 mm,(m=3) 对应 37.97 mm。可以用下面代码把理论位置叠加到数值结果上验证。

orders = -4:4; positions = orders * lambda * f / d; hold on; plot(positions, 0.02*ones(size(positions)), 'ro', 'MarkerSize', 6); legend('FFT 仿真', '光栅理论位置');

结果应该看到数值峰的横坐标和红点几乎重合。峰宽随缝数 (M) 增大而变窄,主极大半宽尺度约为 (\lambda f/(Md)),也就是缝数越多,衍射级越尖锐。如果峰宽明显偏大,先检查M是否正确写进了循环里;如果峰位整体偏移,则需要回看fxx_obs的坐标定义。

4.3 缝宽对缺级的影响

平面光栅衍射中有一个容易忽略的现象:某些级次的主极大正好落在单缝衍射包络的零点上,于是整级消失,称为缺级。判断条件是两个位置相等:

[ \frac{m\lambda f}{d}=\frac{n\lambda f}{b} ]

也就是 (m d = n b)。当 (b/d=1/2) 时,偶数级次会被缺掉。上面的仿真里,(m=\pm2,\pm4) 的峰应该很弱或完全消失,而 (m=\pm1,\pm3) 仍然可见。这个现象用光栅方程本身解释不了,必须回到单缝包络和光栅因子的乘积关系。若把b改成 20 μm,缺级变成 (m=5) 的整数倍,因为 (d/b=2.5);若改成 (b=d),则所有级次都会消失,因为此时光栅退化为无调制的透光区域。

5. 仿真的验证、导出与进阶用法

5.1 与解析解叠加比对,确认坐标缩放系数

做完仿真后不能只凭“形状像”就认定正确。最直接的验证是把解析强度曲线和 FFT 数值结果画在同一个坐标系里。以单缝为例:

beta = pi * a * x_obs / (lambda * f); I_theory = (sin(beta)./beta).^2; I_theory(beta==0) = 1; hold on; plot(x_obs, I_single/max(I_single), 'b'); plot(x_obs, I_theory/max(I_theory), 'r--');

如果两条线几乎重合,说明波长、焦距和坐标轴映射都正确。注意beta==0的位置会出现sin(beta)/beta的除零问题,手工赋 1。双缝也可以用cos(pi*d*x_obs/(lambda*f)).^2 .* I_theory做对照,但别忘了包络调制。这个步骤能很快暴露坐标缩放系数错误,比如漏乘lambda*f,这时数值峰的横坐标会差一万倍左右,一眼就能看出来。

5.2 把运行结果保存为 .mat 和图片,方便课程报告复用

仿真跑一遍很快,但实验结果、参数和图像要一起保存,才能放进报告或后续论文里。推荐在代码末尾加入:

save('diffraction_results.mat', 'x_obs', 'I_single', 'I_double', ... 'I_grating', 'lambda', 'f', 'D', 'N'); imwrite(uint8(255*mat2gray(I_grating)), 'grating_orders.png');

mat2gray会把强度数据归一化到 0 到 1,再乘以 255 转成 8 位灰度图。对一维强度写成 png 其实意义有限,主要便于归档;真正的定量参数还是要存.mat。如果要进一步做二维仿真,把fft改成fft2,用meshgrid生成xy,再通过imagesc(x_obs, y_obs, I)展示二维图样。这里同样要注意ifft2fftshift的配对,二维情况下调用fftshift时默认只对第一维操作,需要写成fftshift(fftshift(E_far,1),2),或者直接用fftshift(E_far,2)指定维度,否则图样中心会上下错位。对于超大规模采样,比如N=8192,建议把E0转成single,内存占用能省一半,FFT 本身也会更快。

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

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

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

立即咨询