单像素成像之傅里叶频谱采集：从原理到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_{x=0}^{N-1} \sum_{y=0}^{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 次投影，重建质量在视觉上几乎没损失。

这个表格是仿真阶段非常实用的参考量。我在实际项目里往往从 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 误差扩散保留的条纹连续性更好，低频噪声也更少。

但要注意，抖动后的图案调制深度会下降，不再等于你设定的 ( b=0.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.5+0.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 我个人的三个小习惯

最后分享几个踩过坑以后养成的习惯。第一，每次改动代码或硬件参数，都在存档文件名里写上采样率、图像尺寸、有效调制深度这三个数，否则三天后回来看数据，完全想不起是哪个配置跑出来的。第二，测量原始数据和重建代码分成两个独立部包，测量数据永久保留，因为你永远不知道后面改算法后会不会重新解算一遍。第三，环境光是不是稳定，直接决定你低频重建的底噪，我每次实验前都会先测一组黑底电平，若波动超过阈值就先不急着做测量，先把光学平台稳定下来再说。

这套傅里叶频谱采集的单像素成像方案，从仿真到硬件都不复杂，但每一步都有值得深挖的细节。如果这篇过程梳理能帮你少走几个弯，那这几千字就没白写。