简介：这份资源面向学习大学物理光学课程、需要借助编程理解衍射与干涉原理的学生和教师，用 MATLAB 复现经典光学实验的计算机模拟过程。文档围绕两条技术路线展开：一是衍射积分方法，通过计算接收屏上的光强分布，模拟单缝衍射、杨氏双缝干涉与黑白光栅衍射；二是傅立叶变换方法，利用屏函数功率谱模拟简单孔径、黑白光栅及正弦光栅的夫琅和费衍射，并给出可直接运行的完整代码与对应光强分布图。资源包共 1 个 doc 文件，约 247KB，内容涵盖理论公式推导、程序清单与结果图示，便于对照理解衍射因子调制、缺级现象等细节。目前已有 793 人学习，适合作为光学实验验证与课程设计的参考材料，帮助读者把抽象公式转化为可视化的衍射图样。
1. 光学衍射 MATLAB 模拟：为什么两套算法能算出同一张衍射图
很多人第一次看到「用 MATLAB 模拟光学衍射」的代码时，会有一个直觉上的困惑：单缝衍射明明可以用菲涅耳-基尔霍夫衍射积分一个点一个点地累加出来，为什么换个孔径，教材就突然改用fft2了？更奇怪的是，两种方法算出来的图样居然对得上。这背后不是巧合，而是傅里叶光学里一个相当硬的结论：夫琅和费衍射场的复振幅分布，本质上就是孔径函数的傅里叶变换。搞懂这条等价关系，你才知道什么时候该老老实实做积分、什么时候可以一把fft2甩过去，也才知道为什么光栅会出现缺级、为什么正弦光栅只有三个斑。这份 MATLAB 资源把两条路线都写成了可运行脚本，正好拿来当解剖样本。
2. 衍射积分法：把连续孔径拆成有限个相干点源
2.1 从惠更斯原理到离散求和
衍射积分法的出发点非常朴素：把孔径看成一大堆彼此相干的点光源，屏幕上某一点的复振幅，就是这些点源传到该点后叠加的结果。对单缝来说，缝宽为a，分成np个点，第k个点源的位置是yp(k)，屏幕上的观察点纵坐标为ys(i)，距离为D。在傍轴近似下，光程差近似为yp*sinφ，其中sinφ = ys/D。
于是每个点源贡献的相位是alpha = 2*pi*yp*sinφ/lam，把所有点源的cos(alpha)和sin(alpha)分别累加，再平方求和，就得到该点的相对光强。资源里的核心循环是这样的：
for i = 1:ny sinphi = ys(i)/D; % 傍轴近似下的衍射角正弦 alpha = 2*pi*yp*sinphi/lam; % 每个点源相对中心的相位 sumcos = sum(cos(alpha)); % 实部累加 sumsin = sum(sin(alpha)); % 虚部累加 B(i,:) = (sumcos^2 + sumsin^2)/np^2; % 归一化光强 end这里有几个参数值得掰开看。np=51是离散化点数，它决定了积分的精度：点太少，求和结果会偏离连续积分，旁瓣位置会漂；点太多，循环变慢但精度提升有限，因为误差主要来自傍轴近似而不是离散化。ny=51是屏幕采样点数，它决定图样的纵向分辨率。ym=3*lam*D/a是屏幕半宽，取了三倍中央亮纹半宽，刚好能把前几级暗纹包进来。/np^2这个归一化因子不是随便写的：np个点源同相时，实部累加为np，平方后是np^2，除以它就把中央最大值归一到 1，方便后面映射到 0–255 的灰度。
提示：
alpha里的2*pi和后面光栅公式里的pi不是笔误。单缝这里用的是相位差2πΔ/λ，而光栅公式里的α、β是半相位差，差一个因子 2，推导时别混。
2.2 杨氏双缝：为什么光强写成 4cos²(φ/2)
双缝的代码没有走「点源累加」那条路，而是直接用了两束相干光的干涉公式。两条缝到屏幕上 P 点的距离分别是r1和r2，相位差φ = 2π(r2-r1)/λ，光强I = 4I0 cos²(φ/2)。这个4来自两束等强度光相干叠加：振幅相加后是单束的 2 倍，光强就是 4 倍。
for i = 1:n r1 = sqrt((ys(i)-a/2).^2 + D^2); % 上缝到观察点距离 r2 = sqrt((ys(i)+a/2).^2 + D^2); % 下缝到观察点距离 phi = 2*pi*(r2-r1)./lam; % 相位差 B(i,:) = sum(4*cos(phi/2).^2); % 干涉光强 end注意这里r1、r2用的是精确的欧氏距离，没有做sinφ≈ys/D的傍轴近似。这是双缝和单缝在实现上的一个微妙差别：双缝的缝间距a=2e-3相对D=1不算太小，用精确距离更稳；而单缝的缝宽更小，傍轴近似足够。sum(4*cos(phi/2).^2)外面套了个sum，其实对单个标量没意义，更像是从向量化写法遗留下来的习惯，不影响结果，但读代码时容易让人愣一下。
2.3 光栅衍射：两个因子的乘积与缺级
黑白光栅的公式是I = I0 (sinα/α)² (sin(Nβ)/sinβ)²，其中α = πa sinφ/λ是单缝衍射因子，β = πd sinφ/λ是缝间干涉因子，a是缝宽，d是光栅常数，N是缝数。资源里取N=2、d=5a，也就是每条缝等宽、缝间距是缝宽的 5 倍。
alpha = pi*a*sinphi/lam; % 单缝衍射因子 beta = pi*d*sinphi/lam; % 缝间干涉因子 B(i,:) = (sin(alpha)./alpha).^2 .* (sin(N*beta)./sin(beta)).^2;这个乘积结构是理解光栅图样的钥匙。干涉因子决定主极大的位置，条件是d sinφ = mλ；衍射因子决定每个主极大的包络强度，条件是a sinφ = kλ时包络为零。当某个干涉主极大的位置恰好落在衍射因子的零点上，这一级主极大就消失了，这就是缺级。资源里d=5a，所以第 5、10、15……级主极大缺级。你在图上看到主极大之间强度被包络压下去，不是画错了，是物理本身如此。
注意：
sin(N*beta)./sin(beta)在beta接近nπ时是 0/0 型，MATLAB 会给出NaN。资源里靠采样点恰好避开这些点才没翻车，但如果你改n或ym，很可能在图上看到一条NaN造成的断线。稳妥写法是加一个极小量eps到分母，或者用sinc类函数重写。
3. 傅里叶变换法：为什么夫琅和费衍射就是功率谱
3.1 从屏函数到频谱：fft2 的三步处理
傅里叶方法的逻辑是：夫琅和费衍射场的复振幅正比于孔径函数的傅里叶变换，强度就是它的模平方，也就是功率谱。所以只要把孔径图像读进来，做二维 FFT，取模，再搬移零频分量，就能得到衍射图样。资源里的处理链条是imread → fft2 → abs → fftshift → imshow。
a = imread('aperture.bmp'); % 读入孔径图像，灰度矩阵 afft = fft2(a); % 二维离散傅里叶变换，得到复数频谱 aabs = abs(afft); % 取模，得到振幅谱 aabss = fftshift(aabs); % 零频搬到中心，符合观察习惯 imshow(aabss, []); % 自动缩放显示fftshift这一步经常被新手忽略，但它是理解「为什么图样看起来是反的」的关键。fft2输出的零频分量在矩阵的四个角上，直接显示会看到中心是暗的、四角是亮的，和实际衍射图样对不上。fftshift交换第一三象限和第二四象限，把零频搬到中心，图样才和实验观察一致。abs取模之后，相位信息丢了，只剩振幅谱，这对看强度分布足够了，但如果你想做相位恢复或全息模拟，就得保留复数频谱。
3.2 孔径形状如何决定衍射图样
资源里用同一套fft2流程跑了圆孔、矩孔、三角孔、单缝、黑白光栅、正弦光栅，甚至字母符号。这正好验证了傅里叶方法的威力：换孔径只需要换输入图像，算法一行不用改。圆孔得到艾里斑，矩孔得到二维 sinc 函数，三角孔得到带六重对称性的图样，黑白光栅得到离散的主极大点阵，正弦光栅只有三个斑。
正弦光栅只有三个斑这件事值得多说一句。正弦光栅的透过率是cos型，傅里叶变换后只有零频和正负一级三个分量，所以衍射图样就是三个亮点。而黑白光栅是方波型透过率，傅里叶展开有无穷多奇次谐波，所以主极大有很多级。资源里先用cos(xs)+1生成正弦光栅图像存成 BMP，再读回来做 FFT，这个「先生成再读入」的流程虽然绕，但恰好说明了傅里叶方法对输入形式不挑：矩阵直接给也行，图像文件也行。
3.3 采样与补零：FFT 模拟里最容易被忽略的坑
用fft2模拟衍射，有一个隐含假设：输入图像的像素尺寸对应孔径的物理尺寸，FFT 输出的频率间隔对应屏幕上的空间频率。如果你只是把一张 100×100 的 BMP 丢进去，得到的图样在物理尺度上是没有标定的，只能看形状不能读数值。资源里的代码没有显式处理这个标定，因为它主要目的是「看花样」，不是「测条纹间距」。
另一个坑是补零。fft2的频率分辨率由输入矩阵尺寸决定，尺寸越小，频谱采样越粗，细条纹可能被漏掉。常见的做法是把孔径矩阵用零填充到更大的尺寸再做 FFT，相当于在频域插值，图样更平滑。但补零不会增加真实信息，只是让显示更好看。如果你要定量比较条纹位置，补零前后主极大位置不变，但旁瓣幅度会有细微差别，这一点在做数值验证时得心里有数。
4. 避坑与排查：那些让衍射图「看起来不对」的细节
4.1 图样上下颠倒或左右镜像
现象：用fft2算出来的衍射图样和教材上的图对比，发现上下或左右反了。原因：fft2的零频在角上，fftshift只搬移了零频位置，但如果你在fftshift之前就做了abs或者用了imshow的默认坐标，显示方向可能和物理方向差一个翻转。解决：统一在fftshift之后显示，并且记住fft2的输入矩阵第一行对应图像顶部，输出频率轴的方向和输入一致，不要额外再flipud。
4.2 光栅图样出现 NaN 断线
现象：黑白光栅的强度曲线上突然缺一段，或者plot出来有一条竖线跳到零。原因：sin(N*beta)./sin(beta)在beta = nπ时分母为零，MATLAB 返回NaN，plot遇到NaN就断线。解决：分母加eps，写成sin(N*beta)./(sin(beta)+eps)，或者用sinc函数重写干涉因子，避免除零。
4.3 单缝衍射旁瓣位置对不上理论值
现象：改小np之后，暗纹位置偏离a sinφ = kλ的理论值。原因：np太小，离散点源不足以近似连续积分，求和结果的零点位置会漂。解决：把np提到 101 或 201，同时保持ny足够大以分辨暗纹；如果还不对，检查ym是否覆盖了足够的衍射角范围，ym太小会把高级暗纹截掉。
4.4 正弦光栅出现多余衍射斑
现象：理论上正弦光栅只有三个斑，但模拟结果出现了五個甚至更多。原因：生成正弦光栅图像时用了cos(xs)+1再量化为 0–255 灰度，量化误差引入了高次谐波；或者图像保存成 BMP 时经过了有损压缩或伽马校正，透过率不再是纯正弦。解决：用浮点矩阵直接做 FFT，跳过图像读写；如果必须走图像流程，保存为无损格式并检查灰度直方图是否只有少数几个值。
4.5 圆孔衍射图样不是同心圆
现象：圆孔衍射的模拟结果出现方形或菱形伪影。原因：输入圆孔图像的分辨率太低，圆的边缘被像素化成阶梯状，高频分量被引入。解决：提高孔径图像分辨率，或者在生成圆孔矩阵时用抗锯齿的软边缘，而不是硬阈值二值化。
5. 两套方法的取舍：什么时候用积分，什么时候用 FFT
衍射积分法和傅里叶变换法不是互相替代的关系，而是覆盖不同的近似区间。积分法从惠更斯原理出发，物理图像清晰，适合教学演示和近场、菲涅耳区的计算，缺点是每个观察点都要对孔径做一次求和，计算量随采样点数平方增长。傅里叶法只在夫琅和费区成立，也就是观察距离远大于孔径尺寸平方除以波长的情况，但它把整个场的计算压缩成一次 FFT，复杂度从 O(N²) 降到 O(N log N)，而且换孔径几乎零成本。
| 维度 | 衍射积分法 | 傅里叶变换法 |
|---|---|---|
| 适用区域 | 近场、菲涅耳区、夫琅和费区 | 仅夫琅和费区 |
| 物理直观性 | 强，逐点累加 | 弱，依赖傅里叶光学结论 |
| 换孔径成本 | 需重写积分表达式 | 换输入矩阵即可 |
| 计算复杂度 | O(N²) 量级 | O(N log N) |
| 相位信息 | 保留 | abs后丢失，需保留复数频谱 |
我一般会这么理解：积分法是「从物理原理硬算」，FFT 法是「用数学定理抄近路」。抄近路的前提是你得确认自己在夫琅和费区，否则近场图样用 FFT 算出来是错的。资源里两种方法都给了，正好可以拿同一个孔径分别跑一遍，对比夫琅和费区的结果是否一致，这是验证自己有没有理解等价关系的最好方式。
还有一个容易被忽略的取舍是相位。积分法天然保留复振幅，你可以直接看相位分布；FFT 法如果只取abs，相位就丢了。做全息、相位恢复、涡旋光这类工作时，必须保留fft2的复数输出，不能只存模。资源里的代码为了显示方便取了模，但你在它的基础上改的时候，记得把afft留着。
从那以后我每次用 FFT 模拟衍射，都会先问自己三件事：观察距离够不够远、输入矩阵的物理尺度标定没有、相位要不要留。这三问不过，图再漂亮也是玄学。希望帮到你。
本文还有配套的精品资源，点击获取
