简介这份资源是面向图像处理与计算成像方向学习者的MATLAB仿真代码聚焦傅里叶叠层成像FPM的图像重建算法实现适合具备一定MATLAB基础、希望理解相位恢复与迭代优化流程的高年级本科生或研究生。压缩包内共1个文件为m脚本类型整体约5KB核心代码围绕低分辨率图像生成、部分傅里叶谱计算、高分辨率图像迭代恢复等环节展开涉及fft2、ifft2等二维傅里叶变换函数以及梯度下降类优化循环并配有图像显示用于验证重建效果。目前已有2600人学习下载说明该方向具有较高的关注度。通过阅读与运行代码读者可以直观理解FPM如何利用多角度照明信息在傅里叶域叠加频谱掌握相位恢复与迭代优化的基本实现思路为后续实验设计、算法改进或相关课题研究提供可复用的仿真框架与调试参考。1. 傅里叶叠层成像到底在解决什么问题从一张模糊图到超分辨重建傅里叶叠层成像Fourier Ptychography简称 FPM这几年在计算成像圈子里热度一直不低核心原因就一个它用一套相对廉价的光学硬件配合 matlab 仿真就能跑出超越物镜数值孔径限制的重建分辨率。你手上如果只有一张低倍物镜拍出来的模糊图传统思路基本没救但 FPM 的做法是换一种采集策略——通过改变照明角度让不同方向的高频信息依次进入系统的通带再用相位恢复算法把这些信息拼回去。这套资源就是围绕这个思路做的 matlab 仿真实现适合做计算成像、相位恢复、超分辨方向的研究生和工程师拿来复现和改参数。它不依赖真实光路纯数值仿真就能验证算法链路对没有光学平台的人来说是刚需。下面我从采集模型讲到代码落地再到踩坑排查把这条链路拆开。2. 傅里叶叠层成像的数学骨架从照明角度到频谱拼接2.1 前向成像模型LED 阵列怎么把高频信息搬进来FPM 的物理基础其实不复杂但第一次接触容易被频谱域的坐标变换绕晕。常见做法是这样用一个 LED 阵列放在样品下方每个 LED 对应一个特定的照明波矢。当第 n 个 LED 点亮时样品出射的复振幅场相当于被一个平面波调制频谱整体平移。数学上写成O_n(k) O(k - k_n)其中O(k)是样品本身的频谱k_n是第 n 个 LED 对应的平移波矢。这个平移量由 LED 位置和样品到 LED 阵列的距离决定。平移之后物镜的有限通带由数值孔径 NA 决定只截取一部分频谱再经过低通滤波和降采样得到一张低分辨率强度图。关键点在于单张图只包含一个子孔径的频谱信息但多角度照明下不同k_n对应的子孔径在频谱域是错开的拼起来就能覆盖更大的频谱范围。我一般会先把这个前向模型在 matlab 里写成函数输入是样品复振幅、LED 位置矩阵、物镜参数输出是一组低分辨率强度图。这一步写对了后面重建才有意义。参数上最容易搞错的是k_n的归一化方式——有人用物理波矢有人用像素坐标下的平移量两者差一个2*pi/lambda的因子混用会导致频谱拼接位置整体偏移重建结果直接糊掉。2.2 相位恢复与频谱拼接为什么需要迭代而不是直接反解前向模型是线性的但问题在于探测器只记录强度相位信息丢了。更麻烦的是每张低分辨率图对应的是频谱域一个子区域的降采样观测直接反解会面临严重的病态性。所以 FPM 的标准解法是交替投影式的迭代先给样品频谱一个初始猜测通常用某一张图的傅里叶变换铺开然后对每个照明角度依次做「截取子频谱 → 逆变换到空域 → 用实测强度替换振幅、保留相位 → 再变换回频域 → 更新对应子频谱区域」。这个循环跑几轮到几十轮频谱的高频部分会逐渐被填上。这里有个反直觉的点迭代次数不是越多越好。我实测下来超过某个轮数后高频区域会开始出现振铃伪影尤其是当 LED 位置标定有微小误差时误差会在迭代中被放大。常见做法是设一个残差阈值残差下降变缓就停而不是死磕固定轮数。另外初始猜测的选择对收敛速度影响很大——用中心 LED垂直照明那张图做初始频谱通常比用边缘 LED 更稳因为中心照明对应的子孔径在频谱原点附近低频信息最完整。2.3 仿真参数表这些数值决定了你能不能复现下面这张表是我在仿真里常用的参数集不同配置下重建质量差异很明显建议先按这套跑通再改参数典型值说明照明波长532 nm绿光常见 LED 中心波长物镜 NA0.1低倍物镜留出超分辨空间LED 阵列间距4 mm决定频谱平移步长LED 到样品距离80 mm影响波矢大小相机像素尺寸1.67 μm决定空域采样率图像尺寸256×256仿真用兼顾速度和细节迭代轮数2050视残差曲线调整频谱重叠率约 30%低于 20% 拼接会出缝重叠率这个参数值得单独说。它由 LED 间距和物镜 NA 共同决定重叠太低时相邻子孔径之间没有共同频段拼接会出现明显接缝重叠太高则采集效率低仿真里表现为需要的 LED 数量暴增。我一般会先算一下理论重叠率再定 LED 阵列规模而不是反过来。3. matlab 仿真代码落地从生成数据到重建出图3.1 生成仿真数据前向模型的代码实现先把前向采集跑通这是后面所有验证的基础。下面这段代码生成一组低分辨率强度图核心是频谱平移加低通滤波% 参数设置 lambda 532e-9; % 波长 NA 0.1; % 物镜数值孔径 pixelSize 1.67e-6; % 相机像素尺寸 N 256; % 图像尺寸 ledPitch 4e-3; % LED 间距 ledDistance 80e-3; % LED 到样品距离 ledNum 5; % 单边 LED 数量共 25 个 % 生成样品复振幅模拟一个分辨率板 [x, y] meshgrid(linspace(-1, 1, N)); sample double(sqrt(x.^2 y.^2) 0.5); sample sample .* exp(1i * 0.3 * sin(10 * x)); % 加一点相位结构 % 物镜通带 fx linspace(-1/(2*pixelSize), 1/(2*pixelSize), N); [FX, FY] meshgrid(fx, fx); pupil double(sqrt(FX.^2 FY.^2) NA/lambda); % 逐 LED 采集 ledRange -(ledNum-1)/2 : (ledNum-1)/2; images zeros(N, N, ledNum^2); idx 1; for i 1:ledNum for j 1:ledNum % 计算波矢平移量 kx 2*pi*ledRange(i)*ledPitch / (lambda*ledDistance); ky 2*pi*ledRange(j)*ledPitch / (lambda*ledDistance); % 频谱平移 spectrum fftshift(fft2(ifftshift(sample))); spectrumShifted circshift(spectrum, round([kx ky] * N * pixelSize / (2*pi))); % 低通滤波 逆变换 spectrumFiltered spectrumShifted .* pupil; field fftshift(ifft2(ifftshift(spectrumFiltered))); images(:,:,idx) abs(field).^2; % 探测器只记录强度 idx idx 1; end end逻辑说明先构造一个带相位结构的样品然后对每个 LED 计算频谱平移量用circshift做整数像素平移仿真里简化处理实际可以用插值再乘物镜通带做低通最后逆变换取强度。参数上ledNum决定采集图像数量ledPitch和ledDistance共同决定平移步长这两个值改了之后重叠率会变重建效果也会跟着变。注意circshift的平移量是取整的仿真里够用但如果要模拟更精细的照明角度得换成频域插值。3.2 重建主循环交替投影的代码结构重建部分的核心是那个交替投影循环下面这段是骨架% 初始化用中心 LED 图像做初始猜测 centerIdx ceil(ledNum^2 / 2); spectrumEst fftshift(fft2(ifftshift(sqrt(images(:,:,centerIdx))))); % 迭代重建 maxIter 30; for iter 1:maxIter for i 1:ledNum for j 1:ledNum % 计算当前 LED 的平移量 kx 2*pi*ledRange(i)*ledPitch / (lambda*ledDistance); ky 2*pi*ledRange(j)*ledPitch / (lambda*ledDistance); shiftPix round([kx ky] * N * pixelSize / (2*pi)); % 截取子频谱 subSpectrum circshift(spectrumEst, -shiftPix) .* pupil; % 逆变换到空域 field fftshift(ifft2(ifftshift(subSpectrum))); % 用实测强度替换振幅保留相位 fieldUpdated sqrt(images(:,:,idx)) .* exp(1i * angle(field)); % 变换回频域更新对应区域 spectrumUpdated fftshift(fft2(ifftshift(fieldUpdated))); spectrumEst spectrumEst circshift((spectrumUpdated - subSpectrum) .* pupil, shiftPix); end end % 计算残差可选提前退出 residual computeResidual(spectrumEst, images, ...); if residual 1e-4 break; end end % 最终重建结果 recon fftshift(ifft2(ifftshift(spectrumEst))); reconAmplitude abs(recon); reconPhase angle(recon);这段代码的关键在更新那一步spectrumEst的更新量是(spectrumUpdated - subSpectrum) .* pupil再平移回去。乘pupil是为了只更新通带内的区域避免把噪声引入到未观测的频段。参数上maxIter设 30 轮通常够用但如果你改了 LED 数量或重叠率可能需要调整。残差计算函数computeResidual需要自己实现常见做法是算重建图与实测图的强度差异。3.3 重建结果评估看振幅还是看相位重建出来之后评估不能只看振幅图。FPM 的一个核心优势是能同时恢复相位所以相位图的质量同样重要。我一般会从三个角度检查振幅图的分辨率是否超过物镜衍射极限、相位图是否连续无跳变、频谱拼接区域是否有明显接缝。如果振幅图清晰但相位图有大量噪点通常是迭代轮数不够或者初始猜测太差如果频谱边缘有亮环多半是重叠率不足导致的拼接伪影。提示仿真里可以用一个已知的相位样品比如透镜相位做 ground truth重建后直接算相位误差比肉眼判断靠谱得多。4. 避坑排查傅里叶叠层仿真里最容易翻车的五个地方4.1 频谱平移方向搞反重建图整体偏移现象重建出来的图像位置和原始样品对不上整体平移了一段距离。原因circshift的平移方向和你计算波矢时的符号约定不一致。有人习惯用kx表示向右平移有人用-kx代码里混用就会导致方向反了。解决在生成数据和重建代码里用同一套符号约定最好把平移量计算封装成一个函数两边都调用它避免手写两遍。4.2 重叠率不足频谱拼接出现明显接缝现象重建图的频谱在子孔径交界处有亮线或暗带空域表现为周期性条纹。原因LED 间距太大或物镜 NA 太小导致相邻子孔径在频谱域没有重叠区域。解决先算理论重叠率公式是overlap 1 - ledPitch / (lambda * ledDistance / (2*NA))低于 20% 就要缩小 LED 间距或增大 NA。仿真里可以直接改ledPitch试。4.3 迭代不收敛残差曲线震荡现象残差跑到一半开始上下震荡重建图越来越糊。原因步长太大或者更新时没有乘pupil做约束导致高频区域被过度更新。解决检查更新公式里有没有乘pupil另外可以加一个松弛因子把更新量乘 0.5 到 0.8 再累加收敛会稳很多。4.4 初始猜测选错收敛速度差一倍现象同样的迭代轮数别人重建清晰你的还是糊的。原因初始频谱用了边缘 LED 的图像低频信息不完整。解决固定用中心 LED垂直照明那张图做初始猜测它的子孔径在频谱原点附近低频最全。这个习惯我每次都会强制走一遍。4.5 相位图出现 2π 跳变看起来像噪声现象重建相位图有大片黑白跳变不连续。原因相位恢复本身有 2π 模糊迭代过程中如果相位更新没有做解包裹处理就会累积跳变。解决仿真里可以在最终输出前做一次相位解包裹matlab 有现成的unwrap函数二维情况下用unwrap(unwrap(phase, [], 1), [], 2)基本能解决。5. 进阶技巧用频谱重叠率反推 LED 阵列参数跑通基础仿真之后最有价值的进阶操作是反过来用重叠率约束去设计 LED 阵列。很多人做仿真时先拍脑袋定 LED 数量和间距跑出来效果不好再回头调效率很低。我的习惯是先定物镜 NA 和波长算出单个子孔径的频谱半径然后根据目标超分辨倍数反推需要的频谱覆盖范围再除以单步平移量得到 LED 数量最后用重叠率下限校验间距是否合理。具体做法假设目标是把 NA 从 0.1 等效提升到 0.3频谱覆盖半径要扩大 3 倍。单个子孔径半径是NA/lambda需要的总平移范围是2*(0.3-0.1)/lambda。如果 LED 间距对应的平移步长是deltaK那单边 LED 数量就是(0.3-0.1)/(lambda*deltaK)向上取整。算完数量后用重叠率公式校验deltaK是否让相邻子孔径有至少 20% 重叠不满足就缩小ledPitch或拉大ledDistance。这套反推流程在仿真里特别省时间因为你不用反复试错。下面这个表是我常用的几组配置可以直接抄目标 NALED 单边数量LED 间距重叠率0.234 mm约 35%0.354 mm约 30%0.473.5 mm约 25%0.593 mm约 22%验证方法也简单重建完之后把频谱估计的幅度图拿出来看有效覆盖范围是不是达到了目标 NA 对应的半径。如果边缘区域幅度明显低于中心说明高频没拼上要么重叠率不够要么迭代轮数不足。我一般会在重建循环里每 5 轮存一次频谱幅度图最后拼成动画看高频是怎么一步步填上的这个习惯帮我定位过好几次收敛问题。从那以后我每次做 FPM 仿真都强制先跑一遍参数反推再写代码省下来的调参时间比写代码本身还多。希望帮到你。本文还有配套的精品资源点击获取