简介面向数字全息研究与光学信号处理初学者的MATLAB代码包聚焦从全息图恢复原始光场的三种典型重构路径基于卷积的近场重构、基于Fresnel变换的中距离重构以及角谱法频域重构。资源共3个文件均为.m脚本压缩包仅1KB代码精简便于逐行阅读和修改参数适合本科毕设、课程实验或入门演练。已有1638人学习浏览常用于对照理论公式理解数字全息重构差异。通过运行和分析这些脚本可以直观比较卷积法、Fresnel衍射积分与角谱法在重构距离、计算复杂度和精度上的取舍同时掌握MATLAB中傅里叶变换、卷积核构造及频域相位恢复的基本实现思路三个脚本模块分别对应近场、中距与远场场景可帮助读者快速搭建自己的全息重构实验框架。 第一次把CCD拍到的干涉条纹变成清晰的重建像时我一度怀疑是代码写错了同一张数字全息图用fresnel函数重建和用conv函数重建输出图像的尺寸完全对不上。后来才意识到这不是错误而是数字全息重构里三种主流算法——菲涅尔变换、卷积法CONV、角谱法——在处理“光怎么传播”这个问题时走了完全不同的三条路。这篇文章就把这三条路各自怎么走、MATLAB代码怎么写、参数怎么定、哪些坑容易踩一次性说清楚。1. 先搞清楚从一张干涉图到一幅重建像计算机做了什么1.1 全息重建本质上是模拟“光往回走”数字全息和普通拍照最大的区别在于CCD记录的只是一张干涉条纹图里面同时混着物光、参考光和它们互相干涉产生的项。重建的目标就是从这张二维强度图里把物光的复振幅场U(x, y)数值解出来。所谓“复振幅”同时包含两个信息振幅对应物体的形貌亮度相位对应物体的高度或厚度变化。用MATLAB做重建本质就是用数值方法求解衍射积分。经典的衍射积分长这样U(x, y) ∬ U0(x, y) · h(x - x, y - y) dx dy其中h是自由空间传播的脉冲响应也就是光从物面传播到像面时每个点发出的球面波在目标面形成的小波函数。这个积分在数学上不好直接算于是研究者们从不同角度做了近似分别得到了菲涅尔变换、卷积法和角谱法三种工程上最常用的离散算法。1.2 三种方法其实是同一个公式换了三条路很多初学者会在这一步纠结“到底哪种方法更高级”。我的看法是它们没有绝对的优劣只是适用的实验条件和重建需求不同。菲涅尔变换近似了球面波的二次相位计算量最小适合记录距离较远的宏观全息。卷积法保留脉冲响应原样通过两次傅里叶变换完成卷积输出像元尺寸和CCD一致。角谱法走的是频域传播路线在近场范围内有严格物理意义是显微镜类全息重建的首选。后面每一章我会给出一段能直接放到MATLAB里跑的代码并解释代码里每个关键行在干什么。2. 菲涅尔变换重建一次FFT搞定但像素大小跟着距离走2.1 算法的数学骨架菲涅尔近似的前提是传播距离d足够大满足菲涅尔近似条件。此时脉冲响应可以写成二次相位因子h(x, y) exp(jkd) / (jλd) · exp[jk/(2d) · (x² y²)]把这个代入衍射积分经过一番整理最终计算式可以拆成一个二维傅里叶变换U(x, y) exp(jkd)/(jλd) · exp[jk/(2d)(x² y²)] · FFT{ U0(x, y) · exp[jk/(2d)(x² y²)] }也就是说先把原始全息图乘以一个二次相位因子再做一次FFT最后外面再乘一个二次相位因子一次 FFT 就搞定所以它叫“单次FFT法”。2.2 一个容易搞反的采样间隔变化用菲涅尔法重建完你一定会发现一个问题输出的图像尺寸和CCD原图完全不一样。原因是输出平面的采样间隔发生了改变dx λd / (Nx · dx)其中Nx是图像横向像素数dx是CCD像元尺寸λ是波长。这意味着重建距离越远输出像素尺寸越大视场也越大。我举个具体数字CCD像元尺寸3.45μm像素1024波长632.8nm重建距离15cm那么输出像素尺寸就是 632.8e-9 × 0.15 / (1024 × 3.45e-6) ≈ 26.9μm。也就是说输出图像的每个像素代表物体平面26.9μm的范围而CCD原始采样是3.45μm一个像素。看起来像是图像被“放大”了视场但代价是分辨率变粗。如果你在小距离下用菲涅尔重建dx会变得很小输出图只有一小块亮斑物体可能只占几个像素这时候就该考虑用下面两种方法。2.3 MATLAB实现function [Urec, dx_out] fresnel_reconstruct(holo, lambda, d, dx) % 菲涅尔变换数字全息重建 % holo : 去除直流后的全息图二维矩阵double类型 % lambda: 波长单位米 % d : 重建距离单位米 % dx : CCD像元尺寸单位米 % Urec : 重建复振幅 % dx_out: 输出平面像元尺寸单位米 [Ny, Nx] size(holo); k 2 * pi / lambda; % 输入平面坐标以数组中心为零点 x (-Nx/2 : Nx/2-1) * dx; y (-Ny/2 : Ny/2-1) * dx; [X, Y] meshgrid(x, y); % 第一个二次相位因子 FFT F fftshift(fft2(fftshift(holo .* exp(1i * k / (2*d) * (X.^2 Y.^2))))); % 输出平面采样间隔 dx_out lambda * d / (Nx * dx); dy_out lambda * d / (Ny * dx); % 输出平面坐标 xo (-Nx/2 : Nx/2-1) * dx_out; yo (-Ny/2 : Ny/2-1) * dy_out; [Xo, Yo] meshgrid(xo, yo); % 第二个二次相位因子 Urec exp(1i * k * d) / (1i * lambda * d) ... .* exp(1i * k / (2*d) * (Xo.^2 Yo.^2)) .* F; end提示fftshift的两处位置是配套的。我见过有人只在外围写一个fftshift出来的重建像中心错位、四角出现半幅图像多半就是这个原因。实际使用中如果菲涅尔重建出来的图像像元尺寸过大、看起来太粗糙可以先把全息图四周补零到更大的尺寸比如从1024补到4096再做fresnel_reconstruct。补零相当于插值能让输出像面采样更密图像更平滑。3. 卷积法重建CONV像元尺寸保持不变但要小心脉冲响应的FFT3.1 为什么卷积法的输出像素和CCD一样大卷积法的思路很直接不对方程做菲涅尔近似直接把衍射积分看成两个函数的卷积然后用卷积定理来计算U IFFT{ FFT(holo) · FFT(h) }这里h就是自由空间的脉冲响应。由于FFT不改变空间域的采样间隔所以输出平面的像元尺寸仍然等于CCD的像元尺寸dx不会像菲涅尔法那样变大。这一点对后续要做像素级测量、定量相位成像的场景非常重要。代价是什么呢你至少需要两次FFT一次对全息图一次对脉冲响应再加上一次IFFT计算量比菲涅尔法大。对于1024×1024这样的尺寸MATLAB跑起来体感差距不明显但原理上确实多算了一两次变换。3.2 最容易踩的坑脉冲响应要ifftshift再做FFT卷积法里最隐蔽的坑是脉冲响应h的坐标排列问题。FFT默认把数组的第一个元素(1, 1)当作原点但我们在生成h的时候用的是以数组中心为零点的坐标网格如果不做处理直接fft2(h)等于把中心在(N/2, N/2)的二次相位函数硬生生平移到了(1,1)频域会多出一块线性相位误差重建结果会是一堆条纹噪声根本看不出物体。解决办法就一行fft2(ifftshift(h))。这也是我在文章开头提到“同样一张图不同算法结果尺寸不同”之外最容易让新手怀疑人生的第二个问题。3.3 MATLAB实现function [Urec, dx_out] conv_reconstruct(holo, lambda, d, dx) % 卷积法数字全息重建 % 输入输出参数同 fresnel_reconstruct % 注意输出像元尺寸保持为 dx [Ny, Nx] size(holo); k 2 * pi / lambda; % 脉冲响应 h坐标同样以中心为零点 x (-Nx/2 : Nx/2-1) * dx; y (-Ny/2 : Ny/2-1) * dx; [X, Y] meshgrid(x, y); h exp(1i * k * d) / (1i * lambda * d) ... .* exp(1i * k / (2*d) * (X.^2 Y.^2)); % 关键一步把 h 的零频从中心移到(1,1)再做FFT H fft2(ifftshift(h)); % 卷积定理 Urec ifft2(fft2(holo) .* H); dx_out dx; end这段代码看起来简单但我建议你在实际项目中做一次自检把holo换成单位脉冲zeros(Nx,Nx); holo(Nx/2,Nx/2)1重建结果应该是一个和h形状一致的亮斑分布。如果出来是斜条纹或者颜色混乱说明ifftshift的位置有问题。3.4 距离特别近时的混叠问题卷积法也不是万能的。当重建距离d过小时脉冲响应h在空间域振荡剧烈采样不足会导致频谱混叠重建图像上出现高频伪影。一个常见补救方法是对h补零到更大尺寸或者用“带限卷积”算法替代直接FFT。如果你只是做本科课设或者验证实验先把距离调大一点通常就能绕开。4. 角谱法重建真正的“精确传播”但要注意倏逝波4.1 从频域角度看光传播角谱法的出发点完全在频率域把物光波分解成无数不同方向传播的平面波每一个平面波在空间中传播d距离后只改变相位、不改变振幅相位改变量由空间频率决定。于是传递函数写为H(fx, fy) exp(j · 2πd/λ · √(1 - (λfx)² - (λfy)²))其中fx、fy是空间频率。整个重建过程就是U IFFT{ FFT(holo) · H }从数学上看只要距离不是太大角谱法被认为是自由空间传播的精确解尤其适合近场全息和显微全息。它的输出像元尺寸也保持为dx。4.2 角谱法和卷积法的关系角谱法里的传递函数H从连续物理上看就是脉冲响应h的傅里叶变换所以角谱法和卷积法理论上应当等价。但在离散实现里这两者的数值结果会有细微差别原因在于卷积法里的H是对有限尺寸的h做FFT得到的受到采样窗口和边界效应的影响。角谱法直接代入解析表达式频率域抽样更干净。因此严格讲角谱法在近场更“接近理论”而卷积法在边界处理上更灵活。实际效果上距离不太近时两者几乎看不出差别但角谱法在近场不容易出现卷积法那种混叠伪影这是它在微观重建领域更受欢迎的原因。4.3 MATLAB实现function [Urec, dx_out] as_reconstruct(holo, lambda, d, dx) % 角谱法数字全息重建 % 输入输出参数同上 [Ny, Nx] size(holo); % 空间频率坐标注意分母是 N*dx不是 dx fx (-Nx/2 : Nx/2-1) / (Nx * dx); fy (-Ny/2 : Ny/2-1) / (Ny * dx); [FX, FY] meshgrid(fx, fy); % 传波数 k 2 * pi / lambda; % 传递函数 H exp(1i * k * d .* sqrt(1 - (lambda * FX).^2 - (lambda * FY).^2)); % 如果出现倏逝波可强制置零避免数值噪声 % H(abs(1 - (lambda*FX).^2 - (lambda*FY).^2) 0) 0; % 和卷积法一样频域函数需要ifftshift到(1,1)位置 H ifftshift(H); % 频域传播 Urec ifft2(fft2(holo) .* H); dx_out dx; end这里有一个很容易忽略的细节空间频率坐标的分母是Nx * dx不是dx。很多人在MATLAB里写频率坐标时习惯写1/dx那是把整个CCD宽度看成采样周期算出来是完全错误的空间频率范围。正确范围是-1/(2dx) ~ 1/(2dx)对应奈奎斯特频率。4.4 倏逝波要不要管当(λfx)² (λfy)² 1时传递函数里根号下变成负数物理上对应倏逝波——它在传播距离极短远小于波长时存在离开表面迅速衰减。在全息重建距离通常远大于波长的情况下倏逝波成分早就衰减没了程序中不处理MATLAB的复数运算也会自动给出衰减结果不会报错。但有些版本会给出复数警告或者在某些异常参数下产生NaN所以我在代码里保留了置零那一行遇到问题可以打开。5. 同一张全息图三种算法跑出来的结果差在哪5.1 我的一组对比测试我用一张模拟的离轴全息图做了三种重建参数如下波长632.8nmCCD像元3.45μm像素1024×1024重建距离150mm。三种方法重建态的强度图肉眼都能看清物体轮廓但细节差别很明显项目菲涅尔法卷积法角谱法FFT次数1次2次2次输出像元尺寸约26.9μm3.45μm3.45μm输出视场大小约27.5mm3.53mm3.53mm重建距离适用范围中远距离近距离为主近、中距离噪声水平低但有边缘振铃近距离混叠风险最干净从表格能直观看出来菲涅尔法适合“看全貌”卷积法/角谱法适合“看细节”。如果你的物体本身只有2mm大用菲涅尔法重建出来只占视场中间很小一块插值放大后模糊换卷积法或者角谱法它就能铺满整个视场轮廓锐利得多。5.2 重建前的“脏东西”直流项和孪生像标题里只提到三种重构算法但实际操作中没有做预处理就直接把原始全息图丢进函数的几乎都不会有好结果。CCD记录的全息图包含三部分直流项零级光、物光真实像、共轭孪生像。如果你用的是离轴全息频谱上这三项是分开的标准做法是对原始全息图做FFT观察频谱。用鼠标或者代码选取1级谱所在区域做二值掩膜。把1级谱搬到频谱中心反FFT得到纯物光波复振幅。对得到的复振幅再做三种算法的传播重建。如果你的实验是近轴同轴全息频谱三项重叠那就只能先用相移法获取物光复振幅再进行传播重建方法不变。我的建议是在调用上述三个函数之前至少先减去全息图的均值把直流项里最亮的那部分削弱。否则重建像中心会有一大块亮斑把物体细节淹没掉。5.3 重建距离怎么选一个简单的自动对焦函数三种方法都有一个核心参数重建距离d。手动一个个试比较费时间我平时会用梯度能量作为评价函数扫描一段距离function score focus_metric(img) img abs(img); % 取强度 gx diff(img, 1, 2); gy diff(img, 1, 1); score mean(gx(:).^2 gy(:).^2); endd_list linspace(0.10, 0.20, 50); % 根据实验距离设定范围 scores zeros(size(d_list)); for i 1:numel(d_list) U conv_reconstruct(holo_clean, lambda, d_list(i), dx); scores(i) focus_metric(U); end [~, idx] max(scores); best_d d_list(idx); fprintf(最佳重建距离: %.4f m\n, best_d);这个方法的原理很朴素聚焦清楚的图像边缘梯度能量最大离焦图像边缘模糊梯度自然下降。对于相位物体可以改用标准差或者拉普拉斯能量但梯度平方和的普适性最好。5.4 关于代码验证再补一句最后分享一个我在教学时经常让学生做的验证实验用MATLAB自己生成一张模拟全息图也就是先设定一个物体的复振幅模拟物光与参考光干涉得到全息图然后用三种算法去重建。如果重建出来的物像能和原始物体对得上说明你的代码链路从头到尾都是对的。模拟全息图的关键代码片段% 模拟物体一个小矩形 obj zeros(512, 512); obj(216:296, 246:266) 1; % 传播距离生成物光用角谱法当“正演” Uo as_reconstruct(obj, lambda, d, dx); % 参考光平面波 Ur ones(size(obj)); % 干涉记录 holo_sim abs(Uo Ur).^2;把这个holo_sim当作holo_clean输入三种算法重建后都应该恢复出那个小矩形轮廓。这一步通过之后再拿实验采集的真实全息图来测试你就能确认问题到底出在算法代码还是光学实验端。我在实际项目里用这套流程处理过颗粒场全息、微结构相位重建也帮学生调过本科毕业设计里的全息重现程序这三段函数经过简单封装基本可以覆盖日常九成以上的数字全息重建需求。工具只是手段搞清楚每种算法背后的坐标约定、采样间隔变化和适用距离才是让MATLAB代码真正“听你话”的关键。本文还有配套的精品资源点击获取