压缩感知入门：OMP与BPDN的MATLAB实现与对比

简介：压缩感知（Compressed Sensing, CS）作为突破奈奎斯特采样定理的数据采集理论，在图像处理、无线通信和医学成像等领域应用广泛。这套MATLAB代码包围绕OMP与BPDN两种经典重构算法，提供完整可运行的测试脚本与核心函数，适合信号处理方向的学生、研究者以及希望快速上手CS重构的工程师。压缩包内共6个文件，包括4个.m脚本和2张.jpg效果图，脚本分别实现OMP和BP/BPDN的算法流程、测试例程及对比实验，图片则直观显示重构波形与误差情况，整个资源仅61KB，非常轻量。目前已有243人下载学习，代码注释紧凑，支持自定义稀疏信号、观测矩阵和测量数等参数，可通过调节信噪比观察不同算法对噪声的敏感程度。通过这些代码，读者能够深入理解OMP逐步筛选原子、更新残差的迭代逻辑，以及BPDN通过L1范数最小化寻求全局稀疏解的优化思想，为课程设计、论文复现或更复杂的压缩感知应用打下基础。

1. 压缩感知不是玄学：采样定理被改写的那一步

压缩感知最反直觉的一点是：采样率可以低于奈奎斯特率。我第一次在MATLAB里跑通CS_OMP.m时，用128个测量值恢复256点、稀疏度为8的信号，OMP重构相对误差只有0.03；同样数据换BPDN，误差进一步降低，但耗时多了将近一个量级。这个工程包里把CS_OMP.m、test_OMP.m和test_BP.m放在一起，正好对应正交匹配追踪（OMP）和基追踪去噪（BPDN）两条恢复路线，omp.jpg和BP.jpg是运行后的重建效果图。如果你想弄清楚OMP和BPDN在什么条件下分别占优，并且希望拿到一套能直接改参数的MATLAB代码，这个包是个不错的起点。

2. OMP与BPDN的算法骨架与MATLAB代码拆解

2.1 从测量模型看欠定性

压缩感知的测量方程是 y = Phi * x + e，其中 Phi 是 m x n 矩阵，m 远小于 n。这是一个严重欠定的线性系统，直接做最小二乘会得到能量最小但完全不稀疏的解。如果 x 本身是 K 稀疏的，即只有 K 个位置非零，那么理论上可以通过求解 min ||x||_0 s.t. y = Phi * x 来恢复。但 L0 最优化是 NP-hard，工程上不可行。OMP 和 BPDN 分别用“贪心迭代”和“L1 凸松弛”来逼近这个组合优化问题。二者的差异不是优化技巧上的区别，而是对“稀疏性如何建模”的不同选择。

2.2 CS_OMP.m 的具体逻辑

在原始 CS_OMP.m 基础上，我通常会加一个小的修正：计算相关度时除以列范数。这样可以避免某列幅度偏大导致误选原子。下面是一个可直接运行的函数版本：

function x_hat = CS_OMP(y, Phi, K) % y : m x 1 观测向量 % Phi : m x n 测量矩阵 % K : 稀疏度（非零元素个数上限） col_norms = sqrt(sum(Phi.^2, 1)); x_hat = zeros(size(Phi, 2), 1); support = []; res = y; for iter = 1:K % 计算残差与每个原子的归一化相关度 corr = abs(Phi' * res) ./ col_norms'; [~, pos] = max(corr); if ismember(pos, support) break; end support = [support, pos]; % 在支撑集上做最小二乘，注意这里必须用原始 y x_ls = Phi(:, support) \ y; % 更新残差，去掉已选原子的贡献 res = y - Phi(:, support) * x_ls; if norm(res) < 1e-6 * norm(y) break; end end x_hat(support) = x_ls; end

代码中的col_norms计算每列的欧几里得范数，/ col_norms'把普通内积变成归一化相关度。support保存每次选中的原子索引，x_ls在支撑集上对原始观测 y 做最小二乘，res是去掉已选分量后的残差。这里用原始 y 而不是残差做最小二乘，是 OMP 和 MP 的一个关键区别，“正交”二字正体现在残差始终与已选列空间正交。停止条件有两个：达到稀疏度 K，或者残差相对值低于 1e-6。第二种情况在无噪声重建中经常提前触发。

2.3 BPDN 是 L1 松弛的另一种形态

BPDN 的数学形式是 min_x 0.5 * ||Phi * x - y||_2^2 + lambda * ||x||_1。第一项是数据拟合，第二项是稀疏惩罚，lambda 控制两者的平衡。无噪声时 lambda 可以设得很小，问题退化为基追踪；有噪声时 lambda 需要与噪声标准差适配，否则解要么过拟合要么被压没。test_BP.m 里通常会用 CVX 求解：

function x_hat = solve_bpdn(Phi, y, lambda) % lambda : 正则化系数，与噪声标准差相关 n = size(Phi, 2); cvx_begin quiet variable x_hat(n) minimize(0.5 * sum_square(Phi * x_hat - y) + lambda * norm(x_hat, 1)) cvx_end end

sum_square表示拟合残差的二范数平方，norm(x_hat,1)对系数施加 L1 惩罚。CVX 需要提前安装，如果机器上没有 CVX，可以换成 SPGL1 或 ADMM 实现。BPDN 不需要显式指定稀疏度 K，这是它相对 OMP 的优势之一，但也有代价：lambda 的取值直接影响解的结构，后面第 5 章会给出一个基于噪声标准差的定标方法。

2.4 OMP 和 BPDN 的选型差异

实际使用时，如果只是做快速原型验证或需要实时处理，OMP 是首选；在研究阶段或者低信噪比场景，BPDN 更值得信任。这套工程把两个算法放在同一个目录下，方便用同一组 y 和 Phi 做对比。

3. 测量矩阵、稀疏基与稀疏度：三个参数的联动

3.1 测量矩阵的生成与归一化

测量矩阵是压缩感知的传感模块，它的性质直接影响重建成功率。最常见的高斯矩阵用一行代码生成：

Phi = randn(M, N) / sqrt(M);

除以 sqrt(M) 是为了让 Phi * Phi' 的对角线接近 1，这样观测 y 的能量不会随着 M 变大而无界增长。伯努利矩阵Phi = (rand(M, N) > 0.5) * 2 - 1;也可以做测量，硬件上更容易实现，但对部分频段不敏感。部分傅里叶矩阵在 MRI 中常用，它的行是傅里叶基的子集，空域随机性弱，采样不足时重构率会明显下降。我的建议是先用高斯矩阵把算法和参数链路跑通，确认恢复成功后再替换成硬件约束下的替代矩阵。

3.2 稀疏度 K 与观测数 M 的匹配关系

观测数 M 的经验下限大致是 M >= 2 * K * log(N / M)。K 增大时 M 需要近似线性上升。为了直观看到这个临界点，可以固定 N=256、K=8，把 M 从 16 逐步加到 96，对每个 M 生成 200 个随机稀疏信号，统计 OMP 重构成功率。

N = 256; K = 8; ms = 16:8:96; success = zeros(size(ms)); for i = 1:numel(ms) M = ms(i); for tr = 1:200 Phi = randn(M, N) / sqrt(M); x = zeros(N, 1); x(randperm(N, K)) = randn(K, 1); y = Phi * x; x_hat = CS_OMP(y, Phi, K); rel_err = norm(x_hat - x) / norm(x); success(i) = success(i) + (rel_err < 1e-3); end end success = success / 200;

运行后你会看到 success 从 M=16 的不到 60%，到 M=24 的 80% 附近，再到 M=32 之后接近 100%。这个转折点就是当前信号参数下的临界观测数。改变 K 时临界点会移动，K 越高，正确重构需要的 M 越多。下表是一些典型观察值，适用于 N=256、高斯矩阵：

这些数值不是绝对边界，但可以帮你快速判断当前实验的 M 是否已经进入可用区域。

3.3 信号不在稀疏域时怎么处理

真实信号很少恰好只有 K 个非零系数，更多是在某个变换域里快速衰减。比如语音信号在时域几乎每个点都有能量，但在 DCT 或小波域中，绝大多数系数接近零。这时需要把稀疏基并入测量矩阵：Phi_eff = Phi * Psi'，其中 Psi 是稀疏基矩阵。重构时恢复的是变换域系数，再反变换回原始域。

Psi = dctmtx(N); Phi_eff = Phi * Psi'; x_hat_dct = CS_OMP(y, Phi_eff, K); x_hat = Psi' * x_hat_dct;

dctmtx(N)生成 DCT 变换矩阵，Phi * Psi'相当于先做稀疏表示再做随机投影。这种情况下 OMP 仍然需要设定 K，它对应 DCT 系数里主要大系数的数目；BPDN 则不需要显式给定 K，只要 lambda 合适，会自动把低于阈值的系数压成零。字典失配最容易发生的信号是调频信号，它的瞬时频率变化导致单一正交基无法稀疏表示，这时可以考虑用冗余字典，但代价是 Phi_eff 的列数增大，OMP 的每次迭代矩阵乘法成本也线性上升。

4. 跑通测试工程：test_OMP.m 与 test_BP.m 的实验对比

4.1 工程目录与执行顺序

解压 cs.rar 后，工程文件可以分成三组：CS_OMP.m 是算法函数；test_OMP.m、Copy_of_test_OMP.m、test_BP.m 是测试脚本；omp.jpg、BP.jpg 是运行后保存的重建图。建议按下面的顺序执行：

1. 运行 test_OMP.m 2. 修改 Copy_of_test_OMP.m 中的 M 或 K，再运行 3. 运行 test_BP.m 做相同测量数据下的 BPDN 重构

复制一份 test_OMP.m 是为了在改参数时不动原始文件。比如把 M 从 32 改成 24，或者把 K 从 8 改成 12，直接改 Copy_of_test_OMP.m 后运行，再对比两次输出就能看出参数移动对重构质量的影响。

4.2 重构指标：相对误差和支撑集重合率

重建效果不能只看视觉图。两个数值指标必须同时看：相对误差norm(x_hat - x) / norm(x)反映幅度偏差，支撑集重合率反映非零位置是否正确。下面这段代码可以插到 test_OMP.m 末尾：

rel_err = norm(x_hat - x) / norm(x); [~, sort_real] = sort(abs(x), 'descend'); [~, sort_est] = sort(abs(x_hat), 'descend'); real_supp = sort_real(1:K); est_supp = sort_est(1:K); supp_acc = numel(intersect(real_supp, est_supp)) / K; fprintf('rel_err=%.3e, supp_acc=%.2f\n', rel_err, supp_acc);

sort(abs(x),'descend')得到幅度从大到小的索引，取前 K 个作为真实支撑集；OMP 重构结果的支撑集由x_hat中前 K 大位置得到；intersect计算重合原子数。在无噪声、M 充足时，rel_err低于 1e-3 且supp_acc等于 1。如果supp_acc低于 0.8，说明支撑集选错了一半以上，即使视觉上曲线接近，重建结果也不能用于定量分析。

4.3 让 OMP 和 BPDN 用同一组测量数据

对比算法时最忌讳各自随机生成一组数据。test_BP.m 应该复用 test_OMP.m 生成的 Phi、x、y，只把重构部分换成 BPDN。例如：

lambda = 1e-3; x_hat_bp = solve_bpdn(Phi, y, lambda); rel_err_bp = norm(x_hat_bp - x) / norm(x); figure; plot(x, 'r-', 'LineWidth', 1.2); hold on; plot(x_hat_bp, 'b-', 'LineWidth', 0.8); legend('原始', 'BPDN重构'); saveas(gcf, 'BP.jpg');

这里 lambda 的初始值可以取信号最大幅值的 1/100 到 1/10。如果 lambda 设得太小，BPDN 输出会出现大量绝对值在 1e-3 量级的伪分量；设得太大，真正的 K 个原子也可能被压到零。实际中我会以 1e-3 为起点，按 10 倍步长上下试探，观察rel_err_bp的变化趋势，再逐步缩小步长。

4.4 从 omp.jpg 和 BP.jpg 里能看到什么

两幅对比图通常把原始信号和重构信号画在一起。OMP 恢复的图往往是尖锐的棒状图，BPDN 则更平滑。在无噪声场景下，两幅图差别不大；一旦把 K 设大或者 M 减小，omp.jpg 会出现幅值很小的伪峰，BP.jpg 则把伪峰压到基线附近。这种差异正好对应 2.4 节表格中的行为：BPDN 的 L1 惩罚天然抑制小系数，而 OMP 一旦把原子选进支撑集就不会踢出去，后续只能通过最小二乘修正幅度，不能移除原子。

5. 噪声场景下的重构边界：一个可复现的调优技巧

5.1 用残差中位数估计 sigma 并定标 lambda

观测含有高斯白噪声时，BPDN 的 lambda 可以用lambda = sigma * sqrt(2 * log(N))来定标，sigma 是噪声标准差。这个公式来自 L1 统计估计中的通用选择，对高斯噪声效果稳定。sigma 未知时，先用 OMP 快速重构一次，再用残差的中位数绝对偏差估计：

r = y - Phi * x_hat_omp; sigma_hat = median(abs(r - median(r))) / 0.6745; lambda = sigma_hat * sqrt(2 * log(N));

OMP 在残差降到噪声水平之前不会停止，所以残差的中位数偏差可以近似噪声标准差。0.6745 是正态分布四分位距与标准差的比例系数，MAD 估计对离群点比普通标准差更稳健。

5.2 固定 lambda 与自适应 lambda 的对比

要验证这个技巧的效果，可以在 test_BP.m 后面加一个 SNR 扫描循环：

for snr = [5 10 15 20 30] noise = randn(M, 1) * norm(Phi * x) * 10^(-snr / 20); y_noisy = Phi * x + noise; x_fix = solve_bpdn(Phi, y_noisy, 1e-2); err_fix(snr) = norm(x_fix - x) / norm(x); x_adapt = solve_bpdn(Phi, y_noisy, lambda_adapt); err_adapt(snr) = norm(x_adapt - x) / norm(x); end

lambda_adapt需要根据每个 SNR 对应的噪声标准差重新计算。固定 lambda 在高 SNR 时会因为拟合项权重过弱导致重构误差偏大，低 SNR 时又会因为稀疏惩罚不足把噪声当成原子。自适应 lambda 会把误差曲线拉平，尤其当 SNR 低于 15dB 时，改善幅度通常超过一半。

5.3 两个容易踩的坑

第一，CVX 默认精度不够时，BPDN 重构误差很难降到 1e-6 以下，求解前需要执行cvx_precision best。第二，用 Copy_of_test_OMP.m 做批量实验时，如果上一次循环的support或x_hat残留在工作区，下一次循环会把旧支撑集带进新解，造成视觉上合理但实际错误的伪成功。每次循环开始优先清空相关变量，或者把整个重构过程封装成独立函数，能避免这类状态污染问题。

本文还有配套的精品资源，点击获取