单像素成像之傅里叶频谱采集：从原理到Matlab仿真与光路搭建

单像素成像听起来像是个反常识的操作探测器只有一枚要画出成千上万个像素的画面。第一次听说这个思路时我也愣了半天但真正动手做了一遍基于傅里叶频谱采集的单像素成像之后才彻底理解——它其实就是把“拍照”这件事换成一堆投影测量和一次逆变换核心根本不在于探测器而在于你怎么调制光场。做这个项目时我重点解决了两个问题单像素探测器的数据流怎么对应到图像频域以及如何在只采集低频子集的情况下重建出视觉上可用的结果最终在 Matlab 里完整复现了这套流程。下面就把这次实现从原理推导到代码细节、再到真实光路搭建会踩的坑全部整理出来。1. 为什么做单像素成像题目背后的思路线1.1 反向思考单个探测器怎么当相机用传统相机靠的是焦平面上的面阵传感器比如 CCD/CMOS 上有几百万个感光单元每个单元负责一个像素。这个方案当然好但一遇到特殊波段就会很尴尬面阵探测器价格翻好几倍分辨率还上不去。比如红外、太赫兹、弱光探测以及一些需要快速响应的工业检测场景单个光电探测器无论是灵敏度还是响应速度都比面阵器件保守很多成本也低不少。单像素成像的思路就是我不需要同时感光几千上万个点而是用一个高性能单像素探测器配合一台能做快速调制的空间光调制器把二维场景“编码”进一串一维光强信号里再用数字方式把图像重建出来。说白了就是把硬件复用率拉满把成像的一部分工作量转移给算法。这个领域里最典型的三个分支是随机散斑关联成像、哈达玛基单像素成像、傅里叶基单像素成像。三者本质上都是“测量场景和一组已知图案的内积”区别在于图案长得不一样重建时用的数学工具也不一样。这次做的傅里叶频谱采集走的正是第三条路。1.2 为什么偏偏选傅里叶而不是哈达玛先说说容易被第一步就被否掉的方案直接随机散斑。用随机矩阵做测量矩阵是最通用的因为它和压缩感知理论兼容性最好看着也高大上。但实际做起来很痛苦——随机矩阵在硬件上不容易有序显示、单次测量的能量效率不稳定重建算法还得上迭代优化动辄几分钟才能算出一张图。哈达玛基要好一点图案是 0/1 二值矩阵而且数学上正交完备重建就是一次矩阵乘法速度很快。但哈达玛基有个致命问题基图案没有频率分层低频信息和高频信息混在一起你要恢复一张连续光滑的灰度图往往得把绝大部分基全测一遍效率并不高。傅里叶基的优势在于频域结构极其清晰。自然界绝大多数图像的能量都集中在低频区域我可以只测低频那一小块频谱就能重建出视觉上能用的图像这就是标题里“高质量高效”的底气。另一个加分项是傅里叶基图案显示的是连续正弦条纹而条纹图案的相位变化对探测器的光强波动不那么敏感加上三相移提取公式可以自动抵消直流分量实验稳定性会好很多。从实现角度看傅里叶频谱采集还有一个隐藏优势重建过程是纯线性、非迭代的一次二维逆傅里叶变换就出图完全不需要调稀疏度、正则化参数这类对新手极不友好的经验值。2. 傅里叶频谱采集核心原理三步相移提取频域的“实部和虚部”2.1 投影采集本质上是在算内积假设待成像场景的光强分布为 ( I(x, y) )尺寸是 ( N \times N )。我生成一组正弦条纹图案 ( P(x, y) )把它显示在空间光调制器上投射到场景然后用单像素探测器测量总的反射光强。忽略掉光学系统的常数因子后探测器读数可以写作[ S \sum_{x0}^{N-1} \sum_{y0}^{N-1} I(x, y) \cdot P(x, y) ]注意这个式子里没有任何点扩散函数这不是传统成像里的卷积而是在算两个二维离散信号的内积。内积的结果是一个标量一次测量只消耗一枚探测器代价是你要做很多次测量才能覆盖足够的信息。再进一步如果图案取正弦形式[ P(x,y) a b\cos(2\pi \frac{f_x x f_y y}{N} \varphi) ]其中 ( a ) 是直流偏置( b ) 是调制深度( \varphi ) 是相位。那么探测器读出的 ( S ) 中除了直流项之外还包含了场景傅里叶分量在 ( (f_x, f_y) ) 频率位置上的实部和虚部信息。遍历不同的 ( (f_x, f_y) )就相当于在频域里逐点“采样”最终把频域网格填满通过逆傅里叶变换就能恢复图像。2.2 三步相移怎么避开直流干扰理想情况下我想从测量值里直接提取出该频点的余弦分量和正弦分量。但实际测量中图案里的直流偏置 ( a ) 会混进读数直接做一次测量根本拆不开直流和信号分量。这里的关键是相移法。对同一个频率点我依次显示三幅相位分别为 ( \varphi 0 )、( \varphi 2\pi/3 )、( \varphi 4\pi/3 ) 的正弦条纹得到三个测量值 ( S_0, S_1, S_2 )。由三角函数展开可得[ S_\varphi aN^2I_{\text{DC}} b(C\cos\varphi S\sin\varphi) ]其中 ( C \sum I \cos\theta )、( S \sum I \sin\theta )就是该频率点傅里叶系数的实部和虚部。三个测量值联立求解[ C \frac{2S_0 - S_1 - S_2}{3b} ][ S \frac{S_1 - S_2}{\sqrt{3}b} ]而标准二维傅里叶变换用的是复指数 ( e^{-i\theta} \cos\theta - i\sin\theta )所以对应到频域复系数[ F(f_x, f_y) C - iS ]这个 ( F ) 就是频域网格该位置上的复数。把所有位置的 ( F ) 填进二维矩阵再做一次二维逆傅里叶变换图像就回来了。这里我一开始踩过一个坑拿仿真数据直接套 ( C iS ) 结果重建出来上下颠倒、左右翻转后面才意识到推导时必须严格对齐 Matlab 里fft2的符号约定。别小看这个负号实验里一旦搞反后面所有环节都会跟着错。2.3 为什么用三相移而不是四相移常见文献里也有四步相移法相位取 0、π/2、π、3π/2每频点要测 4 次。四步相移的公式很好看但它要的测量次数多实际实验中时间成本直接乘以 1.33。三步相移法每个频点只测 3 次在同样图像尺寸和采样率下速度提升明显。不过也要诚实说一下三步相移对相位的精度要求稍高DMD 显示抖动误差或触发延迟都会让相位不准。其实两者的噪声鲁棒性差异没有想象中大多数成熟的单像素成像系统都在用三相移因为速度和存储都在优势侧。我在实际仿真中对比过加入同样的高斯噪声后三相移和四相移重建的 PSNR 相差不到 0.5dB但时间差是实打实的。3. Matlab仿真流程与完整代码3.1 一次仿真需要哪些环节仿真不能一上来就写代码先把数据流理清楚。完整链路是这样的目标灰度图 → 生成频点坐标 → 对每个被采样的频点生成三幅相移条纹 → 计算内积得到测量值 → 套三相移公式恢复频域复系数 → 填充频域矩阵 → 逆傅里叶变换 → 取实部显示。注意一点仿真里我会真的去做“投影求内积”而不是直接拿fft2(Target)的结果来充数。为什么这样因为真实单像素成像的测量误差来自图案量化、系统噪声、同步抖动。如果仿真直接调fft2那等于是在说“我已知全部频谱”完全没有模拟测量过程代码搬到光路上后参考价值会大打折扣。3.2 完整可运行的 Matlab 复现代码我用 Matlab R2018b 以上版本因为代码里用到了sum(..., all)语法。目标图先用内置的phantom函数生成一个 64×64 的经典头模型方便对比重建效果。先写一个提取单频点系数的函数function F fourier_coefficient(I, N, fx, fy, b) % 三相移法提取单个频点的复系数 % I : 目标图像灰度范围[0,1] % N : 图像宽和高 % fx : 水平方向空间频率 % fy : 垂直方向空间频率 % b : 调制深度实验时常取0.5 [x, y] meshgrid(0:N-1, 0:N-1); theta 2*pi*(fx*x fy*y)/N; P0 0.5 b*cos(theta); P1 0.5 b*cos(theta - 2*pi/3); P2 0.5 b*cos(theta - 4*pi/3); S0 sum(I .* P0, all); S1 sum(I .* P1, all); S2 sum(I .* P2, all); C (2*S0 - S1 - S2) / (3*b); S (S1 - S2) / (b*sqrt(3)); F C - 1i*S; end主程序这一段负责生成频域采样掩膜并重建图像clear; close all; clc; N 64; I phantom(N); I im2double(I); b 0.5; samplingRatio 0.2; % 采样率只采低频圆内频点 % 频域坐标中心在 (0,0) [Ugrid, Vgrid] meshgrid(-N/2:N/2-1, -N/2:N/2-1); dist sqrt(Ugrid.^2 Vgrid.^2); % 按半径从低到高取前 samplingRatio 个频点 sortedDist sort(dist(:)); idx min(round(samplingRatio * numel(sortedDist)) 1, numel(sortedDist)); threshold sortedDist(idx); mask dist threshold; % 逐频点测量 Fmap zeros(N, N); for m 1:N for n 1:N if ~mask(m, n) continue; end fx Ugrid(m, n); fy Vgrid(m, n); Fmap(m, n) fourier_coefficient(I, N, fx, fy, b); end end % 重建中心化 - 逆变换 - 取实部 Irec real(ifft2(ifftshift(Fmap))); % 显示 figure; subplot(1, 2, 1); imshow(I, []); title(Original); subplot(1, 2, 2); imshow(Irec, []); title(sprintf(Recon %.0f%%, samplingRatio*100));这段代码跑起来不超过几秒但信息量很大三层循环代表的是每个频点都要做三次“图案投影”测量这正是单像素成像的真实数据体量。你在仿真中把N改成 128 或 256会明显感受到运行时间随频点数线性上升。3.3 采样率、图像尺寸与重建效果的联动关系傅里叶频谱采集的“高效”是有量化依据的。以 64×64 图像为例总频点数为 4096 个每个频点对应 3 次测量完整采完需要 12288 次投影。如果把采样率降到 20%只采低频区域则只需要 819 个频点、约 2458 次投影重建质量在视觉上几乎没损失。采样率频点数投影图案数10000Hz DMD理论耗时5%2056150.06s10%41012290.12s20%81924580.25s50%204861440.61s100%4096122881.23s这个表格是仿真阶段非常实用的参考量。我在实际项目里往往从 10% 采样率起步先看边缘轮廓和整体灰度分布确认系统正常后再根据需求提升到 30%~50%。不需要一上来就拉满因为高频分量对视觉贡献远低于视频的“分辨率焦虑”。4. 从仿真走向真实光路的工程修改4.1 硬件组合与光路布局仿真跑通只是第一步真正搭光路时你才会意识到算法代码里那 20 行核心逻辑对应的是一整套精密光学平台。单像素成像硬件组合按我的经验可以归成三块第一是调制端主流选 DMD数字微镜器件它能做到上万帧切换速率每个微镜只有 0/1 两个状态第二是光学端需要一个准直光源把图案均匀投射到目标区域再让反射光汇聚到探测器第三是采集端单像素探测器建议用带跨阻放大的硅光电二极管模块输出模拟电压再接入数据采集卡同步信号由 DMD 驱动板提供。光路上最需要留意的是图案和场景的对准。DMD 显示的条纹必须完整覆盖被测区域一旦有遮挡或者投影角度不对重建出来边缘就是黑的或者畸变明显。这个步骤没有任何算法能救只能人工调光路。4.2 DMD 只能翻 0/1灰度条纹怎么显示这是数值仿真和物理实验差别最大的地方。傅里叶基的定义里图案是连续灰度正弦条纹但 DMD 微镜只有“亮”和“暗”没法直接显示灰度。常见做法是抖动二值化。把 0~1 的灰度值和一张抖动阈值矩阵逐像素比较输出 0/1。可以用 Bayer 矩阵也可以用误差扩散法我实际测试下来Floyd-Steinberg 误差扩散保留的条纹连续性更好低频噪声也更少。但要注意抖动后的图案调制深度会下降不再等于你设定的 ( b0.5 )。所以实验里的 ( b ) 不能照搬仿真值应该在系统搭建好后先做一个标定显示全亮图案测光强 ( I_{\text{on}} )显示全暗图案测 ( I_{\text{off}} )再显示余弦条纹用读出的交流幅度反推有效 ( b ) 值。如果不做这一步重建图像的对比度会有系统性偏差。4.3 同步、曝光与噪声抑制的实操经验探测器测的是模拟电压DMD 翻图案是一个一个来的两者必须严格同步。我用的方案是 DMD 驱动板输出帧同步方波同时触发数据采集卡的单次采样。每个图案只采样一次吗不是。我每幅图案会采集多个周期然后取均值具体采样次数要看信号稳定性通常 4~16 次平均。平均次数太多会牺牲速度太少噪声压不下来平衡点需要反复试。这里的噪声来源很杂有 DMD 翻动瞬间的抖动脉冲有环境杂散光有探测器本身的暗电流。为了抑制环境光我一般会在探测器前面加一个窄带滤光片光源用对应波长的激光或 LED。如果不加滤光片50Hz 工频闪烁都可能直接叠进频域直流附近低频重建会被污染得很严重。5. 常见问题与排查实录5.1 重建图像灰蒙蒙、对比度低先检查直流项和调制深度标定。三相移公式虽然能自动消去 ( a )但 ( b ) 如果标错所有频域系数都会缩水最后逆变换出来的就是整体灰蒙蒙的图。你可以在频域矩阵里把直流分量拎出来看如果它明显偏离理论值就差不了是 ( b ) 的问题。另一个常见原因是投影图案的调制范围没有铺满 DMD 的灰度映射。比如你代码里写 0.50.5cos但驱动板把 0~1 线性映射到 0~255某些显卡驱动还会做 Gamma 校正输出灰度和光强就不是线性关系了。建议在实验前专门用光电探测器测一遍灰阶响应曲线建立非线性的查找表反向校正。5.2 重建图像有条纹、断层、伪影这类问题八成出在同步和相位误差上。三步相移要求每幅图案的显示时间完全一致触发信号抖动大测量读数之间就会错位归一性就被破坏。我排查时先做“全频谱重建”就是采样率调到 100%如果全采样后图像还有奇怪条纹说明不是采样策略问题而是时序问题。还有一种容易被忽略的是解码时序和图案编号对不上。DMD 一堆图案在板上自动播放采集卡的触发次数一旦错一拍后面所有数据都会对到错误的图案上重建结果会像车祸现场。遇到这种情况我的做法是在开场放一幅独一无二的全亮或者全暗图案做同步标记软件解析时也同步对齐。5.3 分辨率上不去边缘模糊某一批频点没有覆盖到高频边缘自然发糊。解决方向不是盲目提高采样率而是调整采样策略。单纯按圆形低频掩膜采能保证主体轮廓但边缘锐利度不够实际更稳的是“低频全采 高频随机稀疏采样”在某些压缩感知框架里这种策略能显著提高主观清晰度。6. 下一步拓展与工程习惯6.1 从仿真到物理实验的衔接清单如果你准备照这套流程做真实实验我的建议顺序是先在 Matlab 里用模拟测量数据把上面代码全部跑通再做小尺寸 32×32 的硬件验证最后再上 128×128 或更高分辨率。小尺寸的好处是图案数量少、单次测量时间短排查问题快错误定位也直观。做完基础版后可以加的东西太多了RGB 三通道分时采集实现彩色单像素成像、利用压缩感知做稀疏采样、或者把频域掩膜变成自适应策略。我个人的感觉是第一次做千万别贪多先把黑白 64×64 全流程走通硬件同步和标定建立了信心后面的扩展都是顺手的事。6.2 我个人的三个小习惯最后分享几个踩过坑以后养成的习惯。第一每次改动代码或硬件参数都在存档文件名里写上采样率、图像尺寸、有效调制深度这三个数否则三天后回来看数据完全想不起是哪个配置跑出来的。第二测量原始数据和重建代码分成两个独立部包测量数据永久保留因为你永远不知道后面改算法后会不会重新解算一遍。第三环境光是不是稳定直接决定你低频重建的底噪我每次实验前都会先测一组黑底电平若波动超过阈值就先不急着做测量先把光学平台稳定下来再说。这套傅里叶频谱采集的单像素成像方案从仿真到硬件都不复杂但每一步都有值得深挖的细节。如果这篇过程梳理能帮你少走几个弯那这几千字就没白写。