衍射计算与数字全息MATLAB源码解析：从FFT到全息图重建

简介本资源是《衍射计算及数字全息》教材附录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-FFT1次中等及以上距离计算快、内存占用低适合实时性要求高的情况输出平面的采样间隔随距离变化需要对坐标轴做换算角谱法2次任意距离空域和频域采样间隔一致不会改变像元尺寸适合近距离传播当空间频率超出传播波范围时需要处理倏逝波与截止条件卷积D-FFT2次中短距离容易理解算是对冲孔径卷积的一种直观表达结果尺寸需要裁切运算量比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做完频移。我自己的习惯是统一采用“空域坐标原点在矩阵中心使用fftshiftifftshift配对”的约定这样不容易乱。第一次跑通后建议做一件事把圆形孔径的半径、传播距离分别改小和改大观察衍射图样的变化。例如孔径变小衍射条纹会变得稀疏且范围更大距离增大图样整体会变宽。这种直观感受比背一百遍公式都有用。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里的基础函数抽出来加入自己的参数规范和注释再慢慢扩展新的算法比如多重距离重建、相位解包裹、自动聚焦等。等到手里积累了一套稳定的程序库再回头看最开始那些乱成一团的报错就会觉得每一步都值得。本文还有配套的精品资源点击获取