简介面向光学工程与信息光学方向的学习者这套基于MATLAB的全息成像仿真代码以菲涅尔衍射理论为核心完整演示了全息记录与数值再现两阶段过程。压缩包内共有三个脚本大小仅2KB文件极度精简三个脚本分别围绕光源与物体模型构建、干涉图样计算、频域重建等关键环节展开其中涉及fft2、ifft2、fftshift、meshgrid等常用函数便于理解全息图在频域中的处理方法。已有三百八十五人学习下载。运行代码后可以直观观察干涉条纹如何同时保存物体光波的相位与振幅信息进而通过菲涅尔衍射与逆傅里叶变换在指定距离上恢复原始像同时模块化的代码结构方便修改波长、距离、采样点数等参数既能用于课程实验演示也能作为毕业设计或科研入门时的可扩展模板。1. 全息成像仿真为什么绕不开菲涅尔衍射如果你只是把一张图片叠上随机相位然后做一次傅里叶变换那叫“看起来像全息图”不叫全息成像仿真。真正的全息仿真核心是把光的衍射传播过程用数值方法复现出来——而菲涅尔衍射就是其中适用面最广、实现成本最低的传播模型。CX11_2 这类编号经常出现在课程设计或科研预研里本质就是“在计算机里把物光波传播到全息面再把这个波前编码成全息图最后用数值重建把物体捞回来”。整条链路里菲涅尔衍射既是正演的传播算子也是重建的逆运算基础。这件事对做过光学实验的人有吸引力是因为它把实验室里的隔振平台、激光器和干板换成了内存里的复数矩阵而对纯软件背景的人来说它的门槛在于“波长、距离、像素大小”这些参数怎么换算成离散网格。这篇内容按理论、代码、调参、验证的顺序推进目标是你拿到任意一张灰度图能在 MATLAB 里跑出一个可量化的全息重建结果。2. 菲涅尔衍射的适用范围与离散化前提2.1 从 Rayleigh-Sommerfeld 到菲涅尔近似的两个条件严格来说衍射传播可以用 Rayleigh-Sommerfeld 积分描述它对距离和孔径没有额外限制但数值实现里需要对每个物点做卷积复杂度高。工程上更常用 Fresnel 近似它把球面波的二次项展开并保留到一次项得到U(x,y,z) exp(jkz)/(jλz) · ∬ U₀(x₀,y₀) exp{ jk[(x-x₀)²(y-y₀)²]/(2z) } dx₀dy₀这个公式成立的前提是菲涅尔数满足 F a²/(λz) 不能太大同时观察区域离光轴的距离要远小于 z。具体工程判据是 z³ (π/(4λ)) · [(x-x₀)²(y-y₀)²]²_max。在仿真里最常见的问题是取了太小的传播距离比如 0.1 米配 532nm 波长和 2 厘米孔径展开误差会以相位噪声的形式混进重建结果。2.2 离散网格与采样率的三条经验法则进入离散实现前必须把连续坐标映射到像素网格。物面和全息面都用 (N,N) 的矩阵表示物理尺寸 L 对应采样间隔 Δx L/N。菲涅尔衍射离散化后输出面的采样间隔变为 Δu λz/(NΔx)这跟 FFT 自带的分辨率约束有关——你不能同时任意指定输入和输出的物理尺寸输出尺寸一旦确定采样间隔就被波长和距离锁死。做仿真时我一般遵守三条法则一是 N 取 256 到 2048 之间太小的话边缘纹波把细节吃掉太大则内存翻倍而精度受 FFT 舍入限制二是传播距离 z 要满足 z NΔx²/λ否则输出面采样间隔比输入还小重建图会发生混叠三是孔径区域或物体尺寸要小于阵面的 1/4给衍射场留出扩散空间不然边界截断会在重建图里形成亮环。2.3 为什么优先选菲涅尔而不是角谱法很多人一上来就选角谱法因为它只有一个 FFT 和一次逆 FFT看起来更“对称”。但角谱法隐含周期性边界条件物体如果不在阵面中心或边缘强度没有衰减到零伪影非常明显。而菲涅尔单次 FFT 法S-FFT天然带一个二次相位因子等效于给源场加了缓变包络对边缘不那么敏感。另一个不同是角谱法要求传播距离小、近场高频成分丰富时才好用而菲涅尔法在中等距离和远场表现更稳定。全息成像的仿真距离通常在几十厘米到米级这个区间正是菲涅尔的主场。3. 用 S-FFT 实现正向衍射传播的 MATLAB 基础代码3.1 最小可运行的菲涅尔传播函数以下代码是 S-FFT 法的最简实现适合先跑通流程再按需改距离、波长和分辨率。function Uout fresnel_sfft(Uin, lambda, z, delta_in) % S-FFT: 单次FFT实现菲涅尔衍射传播 % Uin : 输入复振幅NxN 矩阵单位 m % lambda : 波长单位 m % z : 传播距离单位 m % delta_in: 输入面采样间隔单位 m [N, ~] size(Uin); k 2 * pi / lambda; delta_out lambda * z / (N * delta_in); % 输出面采样间隔 % 频域坐标注意fftshift保证正确偏移 fx (-N/2 : N/2-1) / (N * delta_in); [FX, FY] meshgrid(fx, fx); % 传递函数二次相位近似 H exp(1i * k * z) / (1i * lambda * z) ... .* exp(-1i * pi * lambda * z * (FX.^2 FY.^2)); % 先补零到2N减少环绕误差可选 Uout ifft2(fft2(Uin) .* H) * (N * delta_in)^2; end这段代码的关键不在 FFT 本身而在坐标构造。fx用(-N/2:N/2-1)/(N*delta_in)而不是linspace(-1/(2*delta_in), 1/(2*delta_in), N)两者数值差一点点但在和二次相位因子相乘后会引入一个线性相位斜坡导致重建图位置偏移几个像素。这里用meshgrid生成二维频域网格配合fft2以后的默认零频在左上角必须保证FX和FY是从负到正的对称序列。delta_out就是输出分辨率同一张图传播距离越远重建像就越小这是菲涅尔衍射“衍射扩展”的直接体现。3.2 参数范围与物理单位校验运行前先做一个量纲自检假设lambda532e-9z0.5N512delta_in20e-6计算delta_out 532e-9 * 0.5 / (512 * 20e-6) 2.6e-5米也就是 26 微米。物体如果是 256 像素的小字输出面它只占约 134 像素因为衍射扩展让它变小了。反过来说如果你想把重建像放大应该减小delta_in或增大z而不是在重建后做插值。判断传播是否处于菲涅尔区用前面的判据算一下z^3和二次项极值之比。对于 1 厘米孔径、532nm 波长z 只需要大于 0.23 米所以 0.5 米距离完全在适用区间内。如果有报错说矩阵维度不匹配检查Uin是否为方阵——size(Uin,1)和size(Uin,2)不一致时meshgrid(fx, fx)会生成错误尺寸的 H。3.3 常见误用把 S-FFT 结果直接当全息图S-FFT 得到的是传播后的复振幅它有振幅和相位两个分量但物理上能被记录下来的只有强度振幅平方或者干涉以后的条纹。很多第一次做全息仿真的同学直接把abs(Uout).^2当成全息图再对它做一次菲涅尔重建结果重建出来的物体被零级和共轭像夹在中间看起来一团糊。正确做法是引入参考光让物光与参考光干涉产生余弦条纹再对强度做重建或数值解调。下一章直接进入这一步用 GS 算法生成纯相位全息图绕开干涉条纹压缩比的问题。4. 从菲涅尔正演到 GS 相位恢复全息图生成与重建闭环4.1 纯相位全息图的生成思路与衍射极限约束可商用的空间光调制器只能调制相位所以计算全息里最常用的是相位型全息图。生成方法有很多Gerchberg-SaxtonGS算法是最容易理解的一种它利用“输入面强度已知”“目标面强度已知”“二者之间传播满足菲涅尔衍射”这三个约束来回迭代。衍射极限在这里的体现是目标面采样点之间的最小间距不能小于λz/(NΔx)如果你想恢复的高频细节比这个极限还细GS 算法再怎么迭代也出不来它只会把能量摊到邻近像素上形成模糊。function phase gs_fresnel(target, lambda, z, delta_in, iter) % GS算法生成菲涅尔域纯相位全息图 % target : 目标强度图灰度范围0~1NxN % lambda, z, delta_in 与传播函数定义一致 % iter : 迭代次数一般20~50次 [N, ~] size(target); A_target sqrt(max(target, 1e-10)); % 目标振幅避免0 % 随机初始相位 phase 2 * pi * rand(N); U A_target .* exp(1i * phase); for ii 1:iter % 正向传播到全息面 U_h fresnel_sfft(U, lambda, z, delta_in); % 全息面约束振幅置为1保留相位 U_h_phase angle(U_h); U_h exp(1i * U_h_phase); % 逆向传播回物面 U_back fresnel_sfft(U_h, lambda, -z, delta_out_of(U_h, lambda, z, delta_in)); % 物面约束强度替换为目标相位保留 phase_back angle(U_back); U A_target .* exp(1i * phase_back); end phase angle(fresnel_sfft(U, lambda, z, delta_in)); endGS 的收敛问题值得注意当目标图里有大面积暗区时暗区像素的振幅接近于零相位更新在那里失去意义容易形成随机相位噪声重建图的背景会出现颗粒感。一个常用的改良是每几次迭代把暗区相位平滑一下或者把振幅下限抬高到 0.01 而不是 0。此外上述代码把反向传播写成fresnel_sfft(U_h, lambda, -z, delta_out)因为 S-FFT 公式在 z 取负值时正好对应菲涅尔逆变换无需另写逆函数。delta_out要从正向传播的结果中拿不能用输入的delta_in否则重建位置和比例都对不上。4.2 查表法生成菲涅尔相位透镜并叠加图像GS 的缺点是需要几十次迭代实时或大批量仿真时开销大。另一种实践思路是“查表法”把菲涅尔波带片的相位提前算好每个像素的调制量只跟它到中心的距离有关然后用查表把图像灰度映射为波带片的局部偏移。实现上就是先把r² (x-xc)²(y-yc)²做二维网格再算phase_lens mod(k * r²/(2z), 2*pi)最后生成hologram exp(1i*(phase_lens alpha * target))其中alpha控制图像的调制深度。这种方法的优点是对图像内容不敏感帧率可以做到实时缺点是重建像边缘会有轻微的波带片残留纹路。4.3 重建端到端验证全息图到重建像的最小脚本无论用哪种方法生成全息图重建都是同一套流程对全息图做菲涅尔传播到某个距离 z_rec取振幅平方。为了更接近物理实验通常还要在全息图外侧补一圈零模拟有限孔径的衍射。下面脚本包含了从生成到重建的完整闭环% 参数设置 lambda 632.8e-9; % 氦氖激光波长 632.8nm z 0.3; % 传播距离 0.3m N 512; delta 15e-6; % 15um像素间距 % 读取或生成目标图 target zeros(N); target(200:312, 80:130) 1; % 简单矩形条用于验证分辨率 % 生成GS相位全息图 hologram_phase gs_fresnel(target, lambda, z, delta, 30); U_holo exp(1i * hologram_phase); % 数值重建 U_rec fresnel_sfft(U_holo, lambda, z, delta); rec_img abs(U_rec).^2; rec_img rec_img / max(rec_img(:)); % 观察重建区域对比矩形条位置 imagesc(rec_img);重建图里矩形条会在阵面中心附近出现旁边伴有一个对称的暗共轭像这是相位型全息图固有的孪生像问题。如果矩形条被零级亮斑盖住可以给全息图乘上一个离焦球面波因子把重建平面在频域上偏移到角落去但这样也会牺牲部分视场。另一个高频出现的错误是重建图倒置——菲涅尔传播在数学上对坐标做了镜像如果你重建出来目标左右颠倒不要怀疑代码这是正常物理现象要纠正的话在目标图输入前做一次上下翻转。5. 全息仿真调试三板斧参数扫描、散斑量化与冗余编码5.1 传播距离错误的最快定位方法如果重建图完全不清晰先不要动迭代次数或相位算法用一维扫描定位正确距离。做法是生成一个点光源作为输入做正向传播后用强度分布判断距离是否处于菲涅尔区。点扩散函数应该是一个平滑的高斯状光斑它的半高全宽理论值约等于λz/L其中 L 是阵面物理尺寸。如果算出来的光斑尺寸和理论值差超过 20%优先检查delta_in和N这两个参数决定了频域最高空间频率进而影响扩散角的准确性。另一个定位技巧是输入一个棋盘格图案重建时观察棋盘格边缘是否出现振铃。如果振铃只在某一侧出现说明采样率不够需要把 N 翻倍而不是调 z。5.2 GS 迭代不收敛的量化判据与暗区处理GS 算法在迭代过程中可以追踪目标面振幅的均方根误差RMSE定义为sqrt(mean((abs(U_back) - A_target).^2))。一个常见的“伪不收敛”现象是 RMSE 前 10 次快速下降之后纹丝不动但图像依然有颗粒感。这时要用散斑对比度来评估视觉质量取重建图强度计算标准差除以均值对比度数值在 0.4 到 0.6 之间属于正常相位型全息重建超过 0.8 说明散斑太强通常需要增加迭代次数或在暗区加入平滑约束。如果 RMSE 在中间某次迭代突然跳高那不是算法问题而是fresnel_sfft里z的符号没变导致正向和逆向传播实际执行了两次正传——这个 bug 在代码里非常隐蔽有经验的工程师也会栽一次。5.3 用冗余编码消除零级像和孪生像的实用技巧相位型全息图的孪生像可以用“离轴参考光”思路削弱在相位分布上叠一个线性相位斜坡等效于让参考光倾斜入射。这个斜坡的斜率不能太大否则超出全息图的采样带宽重建像被截成两半。我通常的做法是加一个直径为 N/3 的圆孔掩膜把相位斜坡限制在孔内这样既能把共轭像推到视场外又不损失主要信息。另一个思路是时间平均生成多张不同随机初始相位的全息图按顺序播放人眼或积分传感器会把散斑平均掉代价是对静止场景才有效。这个技巧配合 DMD 或 SLM 做动态显示时几乎是标配但纯仿真里也可以用来对比单帧和多帧平均后的重建质量。6. 重建像质量的三项硬指标与成像链路整体验证仿真的最后一步不是出图好看而是用指标说明“这套参数下系统能达到什么水平”。我会在重建平面上固定选三个评价量峰值信噪比PSNR、结构相似度SSIM、散斑对比度。PSNR 直接用重建强度图与目标强度图计算SSIM 则看局部纹理保持度这两个指标如果只升不降说明增加迭代次数或提高采样率没有带来冗余过拟合。散斑对比度上面已说过配合频域能量占比会更完整。做链路级验证时把前面所有环节串成一个脚本输入是灰度图输出是重建图的指标数值。分别改变 z、N、lambda看指标的变化趋势。比如 z 翻倍PSNR 通常会先升后降峰值对应的距离就是这套系统的“最佳传播距离”N 增加到 1024 以上PSNR 增长变得平坦说明采样已经达到衍射极限。用这种扫描方式你能在半小时内摸清一个陌生参数组合的脾气也能在给别人讲解时用一条曲线代替“实验效果良好”这句空话。本文还有配套的精品资源点击获取
