简介压缩感知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 endsum_square表示拟合残差的二范数平方norm(x_hat,1)对系数施加 L1 惩罚。CVX 需要提前安装如果机器上没有 CVX可以换成 SPGL1 或 ADMM 实现。BPDN 不需要显式指定稀疏度 K这是它相对 OMP 的优势之一但也有代价lambda 的取值直接影响解的结构后面第 5 章会给出一个基于噪声标准差的定标方法。2.4 OMP 和 BPDN 的选型差异对比维度OMPBPDN迭代方式逐步选原子凸优化一次求解需要 K 值是作为停止条件否由 lambda 控制噪声抑制能力弱残差会把小原子误选强L1 惩罚对噪声更稳健典型耗时毫秒级秒级维护成本自己维护支撑集依赖 CVX/SPGL1实际使用时如果只是做快速原型验证或需要实时处理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 需要近似线性上升。为了直观看到这个临界点可以固定 N256、K8把 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 从 M16 的不到 60%到 M24 的 80% 附近再到 M32 之后接近 100%。这个转折点就是当前信号参数下的临界观测数。改变 K 时临界点会移动K 越高正确重构需要的 M 越多。下表是一些典型观察值适用于 N256、高斯矩阵稀疏度 K建议最小 M低 M 时的典型现象412~16支撑集随机选错824~32高频小原子丢失1648~64伪峰增多、相对误差偏高这些数值不是绝对边界但可以帮你快速判断当前实验的 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); endlambda_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残留在工作区下一次循环会把旧支撑集带进新解造成视觉上合理但实际错误的伪成功。每次循环开始优先清空相关变量或者把整个重构过程封装成独立函数能避免这类状态污染问题。本文还有配套的精品资源点击获取
简介压缩感知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 endsum_square表示拟合残差的二范数平方norm(x_hat,1)对系数施加 L1 惩罚。CVX 需要提前安装如果机器上没有 CVX可以换成 SPGL1 或 ADMM 实现。BPDN 不需要显式指定稀疏度 K这是它相对 OMP 的优势之一但也有代价lambda 的取值直接影响解的结构后面第 5 章会给出一个基于噪声标准差的定标方法。2.4 OMP 和 BPDN 的选型差异对比维度OMPBPDN迭代方式逐步选原子凸优化一次求解需要 K 值是作为停止条件否由 lambda 控制噪声抑制能力弱残差会把小原子误选强L1 惩罚对噪声更稳健典型耗时毫秒级秒级维护成本自己维护支撑集依赖 CVX/SPGL1实际使用时如果只是做快速原型验证或需要实时处理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 需要近似线性上升。为了直观看到这个临界点可以固定 N256、K8把 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 从 M16 的不到 60%到 M24 的 80% 附近再到 M32 之后接近 100%。这个转折点就是当前信号参数下的临界观测数。改变 K 时临界点会移动K 越高正确重构需要的 M 越多。下表是一些典型观察值适用于 N256、高斯矩阵稀疏度 K建议最小 M低 M 时的典型现象412~16支撑集随机选错824~32高频小原子丢失1648~64伪峰增多、相对误差偏高这些数值不是绝对边界但可以帮你快速判断当前实验的 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); endlambda_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残留在工作区下一次循环会把旧支撑集带进新解造成视觉上合理但实际错误的伪成功。每次循环开始优先清空相关变量或者把整个重构过程封装成独立函数能避免这类状态污染问题。本文还有配套的精品资源点击获取
