简介这是一份基于相干衍射成像CDI的MATLAB模拟实现源码包面向光学成像、计算成像以及信号处理方向的研究生、科研人员与工程师尤其适合刚接触相干衍射成像、希望通过Matlab快速建立模拟链路的学习者。资源包内含2个m文件均为MATLAB函数/脚本从波前传播到迭代重建的关键步骤均有体现压缩包整体仅1KB轻量精简便于直接阅读、调试与二次修改。目前已有223人学习/下载对于入门相干衍射成像仿真的参考价值较为明确。通过该源码可了解衍射传播模型的构建方式、迭代重建过程的循环策略与数据流组织思路同时源码结构清晰适合在此基础上发展自定义相位恢复或成像优化算法节省从零搭建仿真环境的时间尤其适合作为课程设计、课题预研或论文复现的起步代码。1. 只拍到强度却要重建复振幅CDI 模拟到底在造什么数据光学成像里最反直觉的一件事是你丢掉了相位还能把物体还原出来。相干衍射成像CDI记录的是远场衍射强度相位信息在探测器上根本没被存下来但通过对物体施加“有限支撑”这个先验再用交替投影迭代居然能同时恢复振幅和相位。模拟源码的价值就是把实验里最难控制的照明分布、噪声水平、支撑域误差全换成已知量让你在 MATLAB 里一遍遍重跑同一条迭代链路观察哪个参数真正影响收敛。这篇博文面向用 MATLAB 做光学仿真、数字全息、相衬成像或刚接触相位恢复的工程师用一个最小 CDI 模拟拆解源码里的关键环节物理模型怎么落成 fft2、支撑域怎么初始化、HIO 迭代为什么在某个参数下停住以及最后怎么确认重建结果是可信的。2. 相干衍射成像模拟的物理模型与 MATLAB 矩阵化2.1 从夫琅禾费衍射到二维 FFT模拟源码的第一行等式CDI 的前向模型很干净物体被相干光照明后出射波在远场传播探测器平面上的复振幅等于物面出射波的傅里叶变换。写成离散形式就是一句fft2这也是几乎所有 MATLAB 模拟源码的第一行核心运算。I(u,v) | FFT{ P(x,y) · O(x,y) } |²其中O(x,y)是物体的复振幅透过率P(x,y)是照明光斑整个乘积就是入射波与物体作用后的出射波。注意这里有一个新手容易忽略、老手也容易踩的点波长、距离、像素物理尺寸这些参数在离散 FFT 模型里不会显式出现因为它们只影响频域坐标的缩放比例不影响重建算法本身。模拟里把波长距离归一化矩阵索引差 1 格就相当于物理空间里差一个固定尺度。MATLAB 里要写对这一步关键不是公式而是象限原点。fft2默认原点在矩阵左上角物理上我们希望波前原点在矩阵中心所以标准写法是输入先做ifftshift输出再做fftshift。这个配对写反一次重建出来的物体就会带一个线性相位斜坡看起来像整体偏转了。wave obj .* illum; % 出射波场 F fftshift(fft2(ifftshift(wave))); % 中心化二维FFT I abs(F).^2; % 探测器强度这里ifftshift把矩阵中心移到(1,1)供fft2处理fftshift再把频域结果恢复到以中心为原点的布局。后续做逆变换时要完全反向配对否则空间域和频域的原点就对不齐自相关支撑也会算歪。2.2 支撑域与过采样比为什么矩阵尺寸不能随便给CDI 能成立不是靠算法多聪明而是靠一个物理先验物体占据了有限大小的区域称为支撑域support。支撑域外的空间里出射波必须严格为零。交替投影算法做的事就是让当前估计反复在两个约束集合之间来回投影傅里叶域要求振幅等于测量的sqrt(I)实空间要求支撑域外为零。两个约束的交集越来越小解才越来越唯一。这就引出一个参数过采样比。它定义为探测器边长与物体支撑宽度的比值一维物体理论下限是 2二维实际经验也要取到 2 到 4。模拟里最常犯的错误是把物体铺满了整个矩阵比如坐落在 256×256 矩阵里的一个 200×200 矩形块。这种配置下缺失相位信息的自由度太多相位恢复算法基本不会收敛误差曲线会从第 10 次迭代开始就平行下滑再跑 2000 轮也没用。实用的生成规则是物体直径控制在矩阵边长的 0.25 到 0.4 倍让物体周围留出至少一倍直径的零背景。后面章节给出的模拟参数会沿用这个经验值。提示过采样比不是越高越好。支撑域只占矩阵的 1/16 时频域采样确实更密但有效信号能量占比太低噪声影响会被放大。0.3 倍边长附近是模拟里比较稳的区间。2.3 把探测器参数映射成 MATLAB 变量一张参数表实验里的探测器参数看起来和矩阵运算隔得很远但模拟源码里每一个都能对应到一个变量或一行处理。下面这张表是我在模拟里最常用的映射关系也方便把实验导出的数据套进同一套脚本。实验中的物理量模拟中的对应MATLAB 实现方式探测器像素数矩阵边长 NN 256;像素尺寸、波长、距离离散 FFT 的坐标缩放归一化不显式建模光子计数不足泊松噪声poissrnd(I / max(I) * count)探测器动态范围强度饱和阈值min(I, max(I) * 1e-2)中心光束阻挡中心小圆掩膜对I中心半径若干像素置 0读出噪声加性高斯噪声I randn(N) * sigma实验里从相机导出的强度图通常是 CSV 或 TIFF 格式CSV 场景下用readmatrix(data.csv)读进来就是矩阵后面的流程和模拟完全一致。唯一要确认的是数值范围和单位模拟代码里一般会先做一次归一化保证强度峰值不是几万量级的大数避免后续poissrnd或sqrt出现数值溢出。3. 用 MATLAB 源码搭一套 CDI 模拟的最小可行流程3.1 生成复振幅物体与照明光斑模拟的第一步是制造一个“已知真值”的复振幅物体。振幅部分可以用图像处理工具箱里的phantom生成一个类似组织结构的灰度图相位部分叠加一个平滑面形让重建任务同时包含振幅对比和相位延迟。% CDI 模拟参数区 N 256; % 探测器像素数矩阵边长 obj_r round(N * 0.18); % 物体半径过采样比约 2.8 rng(0); % 固定随机种子保证可复现 % 生成复振幅物体幅值 相位 amp phantom(N); % Shepp-Logan 模型作幅值 [X, Y] meshgrid(1:N, 1:N); phase 0.4 * peaks(N); % 平滑相位面形峰值约 2.4 rad obj amp .* exp(1i * phase); % 复数透过率 % 照明光斑平面波全1高斯照明时取消注释下一段 illum ones(N); % illum exp(-((X-N/2).^2 (Y-N/2).^2) / (2*(0.5*obj_r)^2)); wave obj .* illum; % 出射波场phantom的输出范围是 0 到 1peaks函数输出大约在 -6 到 6 之间乘 0.4 后相位幅度约 ±2.4 rad既有明显的相位包裹又不至于让相位梯度大得离谱。物体半径取0.18 * N物体外圈有超过自身直径两倍的零背景过采样比满足要求。照明部分默认用平面波。想模拟聚焦照明时把注释掉的高斯光斑换上去高斯宽度用物体半径的一半重建时支撑域会自动适应照明轮廓这也是模拟相对实验的便捷之处。3.2 前向传播与强度记录前向传播只需要一次中心化 FFT然后取模平方。为了接近真实实验加泊松噪声模拟光子计数统计涨落并加一个中心遮挡来模拟光束阻挡。F fftshift(fft2(ifftshift(wave))); % 远场衍射 I abs(F).^2; % 泊松噪声控制峰值光子数影响信噪比 photon_count 1e4; % 峰值光子数 I_norm I / max(I(:)) * photon_count; I_noisy poissrnd(I_norm); I_meas I_noisy / photon_count * max(I(:)); % 还原到原始量级 % 中心光束阻挡beamstop r_bs 3; [xx, yy] meshgrid(1:N, 1:N); mask_bs (xx - N/2).^2 (yy - N/2).^2 r_bs^2; I_meas(mask_bs) 0;poissrnd的输入是期望光子数输出是带泊松涨落的随机计数值。先把强度归一化到峰值photon_count加完噪声再缩放回去是为了让噪声后的强度仍处于原仿真量级既不过度溢出也不丢失动态范围。中心遮挡半径 3 个像素对应实验中阻挡直射光的小圆盘。遮挡区域后续在迭代里要特殊处理不能把它当成真正的零强度测量值。3.3 自相关支撑初始化从强度图猜物体位置相位恢复不能从空支撑出发有个经典做法是通过强度图的自相关来估计支撑位置和尺寸。自相关的支撑宽度大约是物体支撑的两倍所以先对测量振幅做一次逆傅里叶变换取显著区域再缩小一半。% 自相关得到物体支撑的粗估计 ac fftshift(ifft2(ifftshift(sqrt(I_meas)))); ac abs(ac); th 0.25 * max(ac(:)); blob ac th; % 显著区域约等于物体与其翻转的卷积 % 区域质心作为物体中心等效边长折半作为物体尺寸估计 stats regionprops(blob, Centroid, Area); cx round(stats.Centroid(2)); cy round(stats.Centroid(1)); L_ac sqrt(stats.Area); % 自相关等效边长 obj_d max(6, round(L_ac / 2)); % 物体半径估计 % 生成圆形初始支撑 support false(N); [sx, sy] meshgrid(1:N, 1:N); support((sx - cx).^2 (sy - cy).^2 obj_d^2) true;regionprops来自图像处理工具箱如果环境中没有这个工具箱可以用mean和std找质心也可以固定取矩阵中心作为支撑中心——模拟里物体通常就在中心附近。自相关撑起的区域是物体支撑的自卷积所以L_ac / 2是物体尺寸的合理估计。这个初始支撑不必很准后面 HIO 迭代中会逐步修正。3.4 HIO 迭代主循环与 shrinkwrap 支撑更新核心迭代采用混合输入输出HIO算法配合每隔一定轮数更新一次支撑的 shrinkwrap 策略。HIO 在支撑外的处理比误差还原ER多了一个负反馈项能有效摆脱局部极小是 CDI 模拟源码里最常见的主循环。beta 0.8; % HIO 反馈参数 n_iter 500; rng(1); obj_est sqrt(I_meas) .* exp(1i * 2 * pi * rand(N)); obj_est fftshift(ifft2(ifftshift(obj_est))); % 随机相位初始 err_f zeros(n_iter, 1); h fspecial(gaussian, [5 5], 1.0); % shrinkwrap 平滑核 for k 1:n_iter % 傅里叶域保持测量振幅替换相位 F_est fftshift(fft2(obj_est)); F_upd sqrt(I_meas) .* exp(1i * angle(F_est)); g fftshift(ifft2(ifftshift(F_upd))); % 回到实空间 % 实空间支撑内取 g支撑外按 HIO 规则更新 obj_new g .* support ... (obj_est - beta * g) .* (~support); obj_est obj_new; % 傅里叶域相对误差 err_f(k) norm(abs(F_est) - sqrt(I_meas), fro) ... / norm(sqrt(I_meas), fro); % 每 20 轮更新支撑shrinkwrap if mod(k, 20) 0 blurred imfilter(abs(obj_est), h, replicate); support blurred 0.15 * max(blurred(:)); support imdilate(support, strel(disk, 2)); end end这段循环里两个约束交替生效傅里叶域用测量振幅替换估计振幅、保留估计相位实空间则把支撑域内的值直接替换为g支撑域外的值按obj_est - beta * g更新。beta越大支撑外的抑制越强收敛快但容易振荡beta太小则摆脱局部极小的能力弱。0.8 是多数模拟任务里不用怎么调的默认值。支撑更新用的是高斯模糊后取阈值再做半径为 2 的膨胀避免支撑边界收缩得过于激进。初始支撑偏大的情况下这个策略会在迭代中逐渐把支撑拉近到物体真实边界。注意中心遮挡mask_bs区域的强度是人为置零的不是真实测量。严格处理时该区不应参与傅里叶约束更稳的做法是给这块区域权重 0。简单模拟里直接让sqrt(I_meas)为零也能收敛但重建物体中心会有一个轻微暗斑这是 beamstop 伪影而不是算法问题。4. 相位恢复算法的参数设定、收敛判据与抗噪处理4.1 ER、HIO、RAAR 三种更新的差别在实空间投影这一步不同算法的差异只有一两行代码收敛行为却完全不同。误差还原ER是支撑外直接置零单调去逼近一组可行解但很容易在第一个局部极小值附近停住。HIO 引入了反馈项让支撑外的残余误差反向作用到下一次估计上跳出局部极小的能力明显更强。RAAR 则在 HIO 与反射型更新之间做插值对噪声和高饱和度区域更稳健。算法支撑域外更新规则抗局部极小抗噪声典型参数ER直接置 0弱弱无HIOobj - beta * g强中beta 0.7~0.9RAARHIO 与反射的凸组合中强混合系数 0.9 附近MATLAB 里从 HIO 换到 ER 只是把支撑外的一行换成obj_est obj_est .* support但实际使用时不要单跑 ER常见做法是先跑 100 轮 HIO 让支撑收敛再切到 ER 做最后的平滑细化。RAAR 实现更复杂一些等位相恢复遇到明显噪声时再考虑替换入门阶段把 HIO 调好就够用。4.2 beta、迭代轮数与收敛判据的配置beta 不是越大越好也不是越小越稳。从经验看beta 在 0.7 附近时 HIO 的振荡幅度适中降到 0.5 以下每次迭代对支撑外误差的修正太小500 轮下来误差降得又慢又不彻底升到 1.0 以上误差曲线会出现周期性震荡看起来像在两组解之间反复横跳。迭代轮数的判断要结合误差曲线。傅里叶域相对误差err_f的表达式已经在第 3 章代码里定义正常的收敛过程是前 50 轮快速下降然后进入平缓滑行。如果发现曲线走平后还要继续加轮数去“碰运气”正确做法是换一个随机相位初值重启或者先检查支撑是否过小。% 收敛判定最近100轮相对变化小于0.5%就提前停下 if k 100 recent err_f(k-99:k); if (max(recent) - min(recent)) / max(recent) 5e-3 fprintf(在第 %d 轮达到收敛\n, k); break; end end这个判据只适用于模拟里固定真值的场景。真实实验没有真值可对比只能靠傅里叶域误差走平来判断迭代是否进入稳态。如果连走平都看不到优先检查过采样比和支撑初始化而不是调 beta。4.3 噪声、饱和像素与权重掩膜的抗噪策略泊松噪声是 CDI 模拟里最接近实验实际的噪声模型。一个值得记住的结论是在sqrt(I)域做傅里叶约束而不是在I域做本身就接近泊松噪声的极大似然处理这也是主循环里都写sqrt(I_meas)的原因不是图方便而是有统计依据。对于饱和像素和 beamstop 区域处理思路应该彻底反过来完全不约束让相位自由演化。实现方法是用一个权重掩膜控制傅里叶域“哪些像素被强迫等于测量振幅”权重为 0 的位置跳过约束。weight ones(N); weight(mask_bs) 0; % 中心遮挡区不约束 sat_pix I_meas 0.99 * max(I_meas(:)); % 饱和区也不约束 weight(sat_pix) 0; for k 1:n_iter F_est fftshift(fft2(obj_est)); amp_target sqrt(I_meas) .* weight ... abs(F_est) .* (1 - weight); F_upd amp_target .* exp(1i * angle(F_est)); % 后续实空间更新与第3章相同 end这段代码把“振幅替换”改为“振幅混合”权重为 0 的像素保留当前估计振幅而不是强行拉向测量值。加了权重掩膜后beamstop 伪影和饱和点扩散都会明显减弱代价是有效约束像素变少迭代需要更多轮数。5. 重建质量怎么量化误差指标、伪影识别与 shrinkwrap 自适应5.1 模拟场景下可用的三个硬指标模拟的好处是有真值可以直接算误差。三个指标里归一化 RMSE 反映振幅的整体偏差相关系数反映结构相似度傅里叶域误差则是迭代收敛性的直接依据。代码很短mask support (abs(obj) 0.1 * max(abs(obj(:)))); rmse sqrt(mean(abs(obj(mask) - obj_est(mask)).^2 ./ abs(obj(mask)).^2)); cc_amp corr2(abs(obj), abs(obj_est));corr2返回 -1 到 1 的值0.98 以上可以算高质量重建RMSE 需要看具体物体动态范围一般小于 0.1 说明振幅恢复得不错。值得注意的是只用angle(obj_est)直接对比相位没有意义因为整体相位还存在一个任意常数偏移对比前要先把重建相位减去平均值。5.2 从伪影形态反推参数问题重建结果出现两类常见伪影时别急着换算法。第一类是背景散斑噪声表现为支撑外有零散亮点通常是初始支撑过大迭代中 shrinkwrap 阈值太低没有及时收缩支撑。第二类是物体边缘振铃和条纹通常是支撑过小物体真实边界被截断了重建被迫把能量挤进了一个过小的框。一个快速判断技巧是显示重建振幅的动态范围。支撑过大时最大振幅会明显超出真值支撑过小时物体内部会出现波纹状明暗交替且边缘处尤其明显。对应修法是调整 shrinkwrap 的阈值参数而不是重跑整个模拟。5.3 shrinkwrap 阈值与膨胀半径的调节技巧shrinkwrap 实际只有三个可调量高斯平滑核尺寸、阈值、膨胀半径。核太大会把支撑边界抹得模糊核太小则容易把目标内部噪声识别成支撑区域。推荐从[5 5]核、标准差 1.0、阈值 0.15、膨胀半径 2 起步这四个值覆盖了多数模拟场景。噪声明显偏高时把阈值从 0.15 提到 0.2同时把膨胀半径从 2 增加到 4撑起一段缓冲带支撑会更稳。实际跑模拟时我一般固定 beta 为 0.8观察误差曲线在 200 轮后的走势再决定优先调阈值还是调膨胀半径这比盲目增加迭代轮数有效得多。本文还有配套的精品资源点击获取