简介：本资源是《衍射计算及数字全息》教材附录B配套的MATLAB程序源代码集，面向光学工程、物理电子、信息光学等方向的本科生、研究生及科研人员，用于辅助理解光波衍射建模、数字全息图生成与物场重建等核心算法。压缩包共23个文件，含22个功能明确的MATLAB脚本（.m）和1份PDF说明文档，涵盖菲涅尔衍射、傅里叶全息、相位恢复、光栅调制、LJCM系列典型仿真案例等关键模块，代码结构清晰、注释完整，便于逐行调试与原理验证；压缩包仅265KB，轻量易用。已有1295人学习下载，适用于课程实验仿真、毕业设计建模及科研预研验证。读者可直接运行各脚本复现书中图例与数值结果，结合PDF说明深入掌握衍射积分离散化、角谱传播、全息图编码与重构等关键技术环节，显著提升光学数值仿真的实践能力。 做数字全息和衍射计算这行的人，手边大概率都翻过“衍射计算及数字全息”这套资料。尤其是附录B给出的MATLAB程序源代码，真正用起来之后才会明白，它解决的并不只是“算一个衍射图”的问题，而是把整个数字全息链路——从光场传播、全息图生成，到重构像恢复——用FFT和矩阵运算串成了一套可以直接改、直接跑的工具流。这篇文章我会从这套代码的使用者视角出发，拆一下它是怎么设计的、核心函数背后的物理逻辑，以及我实际调试过程中踩过的那些坑。如果你正准备用MATLAB做衍射计算，或者在数字全息实验里拿到了全息图却不知道怎么重建，这篇内容应该能帮你少走不少弯路。
1. 从“附录B”说起：这套MATLAB源代码到底能干什么
1.1 数字全息实验为什么逃不开衍射计算
数字全息和传统全息的本质区别，在于用图像传感器代替了干板记录，用计算机算法代替了光学再现。而不管怎么变，整个系统始终被一件事贯穿：光的传播过程要用衍射理论来描述。物光从物体表面出发，经过一段距离到达记录面，这个过程不是简单的几何投影，而是光场复振幅的衍射变换。要做数字全息仿真，就必须先在MATLAB里把这段传播过程算出来，这也是整套源代码最核心的价值。
提到衍射计算，很多人第一反应是“直接套公式就行”。但实际写代码会发现，连续域的菲涅尔衍射积分公式是一回事，离散域里用FFT实现又是另一回事。离散化带来的抽样条件、补零方式、坐标轴定义、相位因子斜率，每一样都可能让结果面目全非。附录B这套程序的价值，恰恰在于它把这些离散化细节都处理好了，给出了经过验证的基本函数库，使用者不需要每次从零开始推导离散公式。
1.2 附录B代码的定位和适用人群
我理解的附录B，应该是一组以衍射传播函数为基础的MATLAB程序集，覆盖从菲涅尔衍射、角谱衍射到全息图生成、数字再现的完整流程。它不是那种只讲概念的教学代码，而是偏向工程实践的工具包，函数输入输出都比较清晰，参数用物理单位制直接定义，方便和实验系统对照。
适合用这套代码的人，第一类是刚开始接触数字全息的研究生，需要快速上手仿真，验证自己的实验方案；第二类是已经在做实验、但手里只有全息图、需要一套成熟重构算法的工程师；第三类是做光学计算或光束传播仿真的人，想拿现成的衍射传播函数作为中间件，集成到自己的系统里。对这三类用户来说，附录B最大的价值不是让你背公式，而是让你拿到一套能跑、能改、能验证的起点。
2. 核心算法拆解：FFT计算衍射的门道
2.1 抽样条件与离散傅里叶变换的关系
要在MATLAB里实现衍射计算，首先得把连续的复振幅分布离散成矩阵。假设光场在空域里以采样间隔Δx均匀采样，那么在频域里，最高空间频率就是1/(2Δx)，也就是奈奎斯特频率。这个看起来很简单的关系，在衍射计算里会引发一连串连锁反应。
以菲涅尔衍射的S-FFT方法为例，物平面光场U0(x0, y0)经过距离z传播后，观察面光场可以写成：
U(x, y) = exp(ikz) / (iλz) · exp[iπ(x² + y²)/(λz)] · FFT{ U0(x0, y0) · exp[iπ(x0² + y0²)/(λz)] }
这个式子里最关键的是两个二次相位因子。第一个在频域变换前乘上，第二个在变换后乘上。要保证计算结果不混叠，这两个二次相位的局部空间频率不能超过采样率允许的范围。二次相位exp[iπ(x0² + y0²)/(λz)]的局部频率近似为x0/(λz)，最大值的绝对值出现在x0 = NΔx/2的地方，所以要求：
NΔx/(2λz) ≤ 1/(2Δx)
等价于：
z ≥ NΔx²/λ
这个条件如果不满足，FFT出来的衍射场就会混叠，图像里会出现周期性的假条纹和噪点，而且不是调参就能掩盖的。附录B代码里很多函数都内置了判断或提示，就是让使用者不要随便取参数。
2.2 菲涅尔、角谱、卷积三种算法怎么选
衍射计算的主流实现方式主要有三种：单次FFT的菲涅尔近似、角谱法、以及卷积形式（D-FFT或双重FFT）。附录B里通常会把这几类都给出，方便根据不同的实验条件选择。我把它们的核心特点整理成一张表：
| 算法类型 | FFT次数 | 适用距离范围 | 优点 | 主要限制 |
|---|---|---|---|---|
| 菲涅尔S-FFT | 1次 | 中等及以上距离 | 计算快、内存占用低，适合实时性要求高的情况 | 输出平面的采样间隔随距离变化，需要对坐标轴做换算 |
| 角谱法 | 2次 | 任意距离 | 空域和频域采样间隔一致，不会改变像元尺寸，适合近距离传播 | 当空间频率超出传播波范围时，需要处理倏逝波与截止条件 |
| 卷积D-FFT | 2次 | 中短距离 | 容易理解，算是对冲孔径卷积的一种直观表达 | 结果尺寸需要裁切，运算量比S-FFT大一倍 |
实际选型时，如果目标是模拟全息图的记录过程，物面和记录面尺寸相同、采样间隔相同，用角谱法最方便，因为不需要处理输出平面坐标缩放的问题。如果目标是做远场衍射，或者物体的尺寸远小于传播距离，用S-FFT更高效。而卷积法在校正像差和离散系统建模时更自然，但运行时间会明显增加。
这一块我之前有过一次教训：做近场全息模拟时，一开始用S-FFT，距离取了几毫米，结果出来的光场完全是一团模糊的条纹，后来换成角谱法，图像立刻恢复正常。原因就是模拟的传播距离小于z = NΔx²/λ这个下限，S-FFT方法本身就失效了。所以算法不是越复杂越好，关键是和你的物理场景匹配。
2.3 全息图到重构像的闭环如何实现
数字全息的整个流程，可以简化成三步：模拟或获取物光场，叠加参考光形成干涉图并记录强度；然后把这个强度图作为输入，乘上复现参考光；最后用衍射传播算法将光场传播回去，得到物体的重构像。
在附录B的程序结构里，这一闭环通常分得很清楚。比如生成全息图的函数做的是前两步，物光场可以通过散射模型或读取实验数据获得，参考光可以是平面波或球面波；而重建部分则包含再现光场的模拟和衍射逆传播。这个分离设计的优势在于，仿真和实验数据可以共用同一套重建代码，只要把“输入全息图”这一步切换成实验拍摄的强度图就行。
实操时，重建过程往往会遇到“孪生像”问题。因为全息图记录的是强度，实像和虚像同时存在，通常需要在频域里滤波把其中一者分离出来。附录B代码中一般会预留滤波处理的接口，我在使用时习惯加一步频域高通或带通滤波，把+1级和-1级衍射项分开，重构质量能有非常明显的提升。
3. 实操阶段：把附录B跑通并改成自己的工具
3.1 代码文件结构和调用关系
拿到附录B的源代码，第一步不要急着运行，先看一下文件结构。我建议你把它当成一个小型工具包来组织，至少区分出三类文件：衍射传播基础函数、全息图生成与重建函数、测试和演示脚本。基础函数只管输入输出光场，不涉及具体实验参数；生成与重建函数负责拼装物光、参考光和传播过程；演示脚本则把一组可行的参数跑通，并画出仿真结果。
这种分层和MATLAB本身的函数调用机制很契合。比如可以在基础函数里写一个AngularSpectrum(U0, lambda, dx, z)，里面做两次FFT，返回传播后的光场；然后在生成脚本里调用它。如果后期要接实验数据，只要保证函数输入输出格式不变，内部的物理模型可以随时替换，不会动到上层调用逻辑。
我看过不少初学者直接把所有代码塞进一个大脚本，参数全部写死，后面改一个波长要翻遍全文件。这不是附录B代码的设计初衷。建议拿到源代码后，先按我上面说的方式把函数提取出来，建一个自己的本地工具目录，以后做实验直接往里面加脚本就行。
3.2 用一组参数把衍射计算第一次跑通
我建议的第一次运行，不要直接加载全息图实验数据，先用一个简单物体做仿真验证，比如一个矩形孔径或圆形孔径。这样你心里有预期结果，可以迅速判断代码是否正常工作。下面我给出菲涅尔S-FFT方法的一个最小示例，参数设置是按照常见实验条件写的：
% 参数设置 lambda = 632.8e-9; % 氦氖激光波长，单位米 N = 1024; % 采样点数 dx = 20e-6; % 采样间隔，单位米 z = 0.3; % 传播距离，单位米 % 生成坐标系 x = (-N/2 : N/2-1) * dx; [X, Y] = meshgrid(x); L = N * dx; % 构造物平面光场：圆形孔径 cx = 0; cy = 0; radius = 5e-4; U0 = double((X - cx).^2 + (Y - cy).^2 <= radius^2); U0 = U0 .* exp(1i * 0 * X); % 可以叠加初始相位 % 菲涅尔S-FFT衍射 k = 2 * pi / lambda; fx = (-N/2 : N/2-1) / L; [FX, FY] = meshgrid(fx); H = exp(1i * k * z) / (1i * lambda * z) .* exp(1i * pi * (FX.^2 + FY.^2) / lambda * z); % 注意：这个写法是频域直接计算公式，需要根据具体附录B函数来调整 % 观察面坐标 dx_out = lambda * z / L; x_out = (-N/2 : N/2-1) * dx_out; [Xo, Yo] = meshgrid(x_out); U = ifftshift(fft2(fftshift(U0))) .* H; U = ifftshift(U); % 根据库函数风格决定是否需要上面这段示例里，fftshift和ifftshift的使用要特别小心。附录B里每个函数对坐标轴的处理方式可能不完全一样，有的在开头就把原点放到矩阵中心，有的则用fftshift做完频移。我自己的习惯是，统一采用“空域坐标原点在矩阵中心，使用fftshift+ifftshift配对”的约定，这样不容易乱。
第一次跑通后，建议做一件事：把圆形孔径的半径、传播距离分别改小和改大，观察衍射图样的变化。例如孔径变小，衍射条纹会变得稀疏且范围更大；距离增大，图样整体会变宽。这种直观感受，比背一百遍公式都有用。
3.3 针对实验数据修改参数的关键位置
如果要从仿真切换到实验全息图重建，需要改的参数就不是简单的波长和距离了。我觉得至少有三个地方必须仔细核：
第一是像素尺寸。相机的像元尺寸决定了采样间隔Δx，常见比如2.2微米、3.45微米、4.65微米等。这个值必须和实验相机一致，否则重建像的物理尺寸会整体漂移。第二是记录距离。这个距离在实验里不是随便估的，通常要测光路的光程差，或者通过自动聚焦算法来搜索。第三是参考光的参数，比如平面波入射角，或者球面波的焦点位置，这些参数会直接决定重现像中心位置和像差。
附录B的代码里，这些参数往往会集中在脚本头部或者在函数入口处定义，这是设计上很友好的地方。不过我要提醒一点：不要直接沿用示例里的默认值。数值稍有偏差，重构像可能看起来还像那么回事，但当你需要量测物体尺寸或做相位定量分析时，误差就会被放大。我在做显微镜数字全息时，就曾因为像素间距少写了一个数量级，导致重构像整体缩放了十倍，一开始还以为是算法错了。
4. 常见问题排查与避坑记录
4.1 输出图像全黑或全白
遇到这种情况，十有八九是数据范围的问题。FFT计算出来的复振幅动态范围极大，直接取实部或者取模，很可能会超出图像的显示范围。如果直接用imagesc(abs(U))而不加坐标轴归一化，亮区会占据全部显示范围，弱信号全被压到黑色里，看起来就是全黑或全白。
解决办法很简单：显示前先对幅度做归一化，或者取对数压缩动态范围。比如：
amp = abs(U); amp = amp / max(amp(:)); imagesc(x_out, x_out, amp); axis image; colormap gray;另外还要检查是不是相位和幅度没有分开处理。数字全息重构后，我们要看的通常是强度分布IG = |U|²，而不是复振幅的实部。有人直接把real(U)拿来显示，结果正负相消，图形完全乱掉，这也是一个很低级但很常见的错误。
4.2 再现像位置偏移、尺寸不对
重构图像位置偏移，最常见的原因是频域滤波中心和参考光角度设置不准。数字全息图在频域里有零级项和两个共轭像项，如果重建时没有把参考光的入射角补偿掉，像的位置就会偏离视场中心。附录B代码里通常会有angle参数或者滤波窗设置，你要做的就是把频域里那对衍射项的坐标找到，让滤波中心对准它。
尺寸不对则多半是采样间隔换算的问题。S-FFT方法里输出平面的采样间隔是dx_out = λz/(NΔx)，和输入平面的采样间隔不相等。如果你忘了这个换算，仍然用原来的坐标轴去画结果，图像看起来就会被拉长或者压缩。用角谱法时输出间隔等于输入间隔，不需要换算，这也是我为什么推荐近场计算尽量用角谱法的原因。
4.3 参数单位混乱导致结果失真
MATLAB代码里单位问题是一个经典陷阱。波长用微米还是米、距离用毫米还是米、像素尺寸用微米还是米，只要有一个对不上，整个系统的比例就会错。我的习惯是全部统一到国际单位制，波长、间距、距离都用米，最后输出图形时再换算成毫米或微米做标注。
还有一点，角度参数容易被忽略。比如平面参考光的角度如果被写成弧度制的，而代码里意外用度数代入，那全息图频域的两个衍射项会偏离得很厉害，导致重构时根本找不到像。这种问题不报错，只体现在结果的形态上，排查起来非常费时间。建议在脚本头写清楚每个参数的单位，并加注释。
4.4 数据太大跑不动？优化思路
数字全息图像通常是千万像素级别，比如2048×2048甚至更大。直接对全尺寸做FFT，速度会变慢，内存占用也很大。附录B的基础函数一般写得比较朴素，内部用的是MATLAB原生fft2，并没有做太多优化，所以跑大数据时可以考虑几个思路。
一是把重建前做频域裁剪，只保留包含目标像的感兴趣区域，而不是对整个频域平面做逆变换。二是利用阈值或下采样降低数据量，前提是不影响你要观测的空间频率范围。三是把多次重复计算的操作写成MEX函数，或者用GPU加速，MATLAB的gpuArray能直接支持FFT，速度提升明显。不过我个人的体会是，先确认算法和参数没问题，再做性能优化，否则在错误的结果上加速没有意义。
5. 最后再分享一点个人体会
这套附录B的MATLAB源代码，我前前后后改过好几轮，用来做过孔径衍射、会聚球面波模拟，也处理过实验记录的全息图。我的感受是，它最大的价值不在于“开箱即用”，而在于给你提供了一套可以对照的离散化范式和函数边界。你遇到问题翻看它的实现方式，能学到很多“课本不会写、但工程必须处理”的细节，比如坐标轴翻转、相位因子的符号约定、以及为什么有的地方要用fftshift而有的地方不用。
如果你准备长期做数字全息相关研究，我建议在这个基础上构建自己的工具箱，不要每次从零写。把附录B里的基础函数抽出来，加入自己的参数规范和注释，再慢慢扩展新的算法，比如多重距离重建、相位解包裹、自动聚焦等。等到手里积累了一套稳定的程序库，再回头看最开始那些乱成一团的报错，就会觉得每一步都值得。
本文还有配套的精品资源，点击获取
