资讯动态

压缩感知MATLAB实战:CoSaMP、SAMP等五种重构算法代码解析与避坑指南

发布时间:2026/10/26 22:57:36 来源:尧图企业网站定制
简介这份资源是面向信号处理、图像重建与无线通信方向的学生和研究人员整理的压缩感知重构算法MATLAB代码合集旨在解决从少量测量值中恢复稀疏信号时缺少可对比实现的问题。压缩包共5个文件全部为m脚本整体约7KB分别对应CoSaMP、SAMP、SP、ROMP等经典重构算法的具体实现便于逐一阅读、运行与横向比较。已有149人学习关注说明其在相关课程实验与课题入门中具有一定参考价值。读者可借助这些脚本理解各算法在支撑集选择、残差更新与收敛速度上的差异对比不同稀疏度与测量数下的重构表现并在此基础上调整参数或设计新的重构策略为压缩感知相关项目提供可直接复用的代码基础。1. 从 cs-code.zip 说起为什么这五个 MATLAB 文件值得你花一晚上跑一遍如果你正在做压缩感知CS相关的课题或者被导师要求“复现一下 CoSaMP 和 SAMP 的对比实验”那你大概率已经在 GitHub 和 MATLAB File Exchange 之间翻了很久。cs-code.zip这个包不大里面就五个.m文件CS_CoSaMP.m、CoSaMP.m、CS_SAMP.m、CS_SP.m、CS_ROMP.m。它没有花哨的 GUI也没有冗长的说明文档但恰恰是这种“裸代码”最适合拿来拆解算法内核。压缩感知重构算法要解决的问题很明确给定观测向量 ( y \Phi x )其中 ( \Phi ) 是 ( M \times N ) 的测量矩阵且 ( M \ll N )如何从 ( y ) 中恢复出稀疏信号 ( x )。这五个文件覆盖了贪婪类算法中最常被引用的几个变体——CoSaMP、SAMP、SP、ROMP以及一个带CS_前缀的 CoSaMP 实现。适合谁适合已经读过压缩感知入门教材、能看懂伪代码、但自己从零写又总觉得“差一点”的硕士生和初级算法工程师。接下来我不谈理论综述直接按“怎么跑通、怎么改参数、哪里会翻车”的顺序拆。2. 五个 .m 文件的角色分工先搞清楚谁调用谁2.1 文件名里的CS_前缀不是装饰打开压缩包第一眼会注意到命名有点乱既有CoSaMP.m又有CS_CoSaMP.m。常见做法是不带前缀的CoSaMP.m是算法核心函数输入测量向量y、测量矩阵Phi、稀疏度K输出重构信号x_hat和支持集。而CS_CoSaMP.m往往是一个外层包装或测试脚本负责生成随机稀疏信号、构造高斯测量矩阵、调用核心函数、计算重构误差并画图。同理CS_SAMP.m、CS_SP.m、CS_ROMP.m大概率也是“测试调用”的复合脚本。我一般会先在工作区执行whos看变量再对每个文件用edit打开扫一眼函数声明行确认哪些是function开头、哪些是脚本。这一步花两分钟能避免后面调用时出现“未定义函数”或“输出参数太多”的低级报错。2.2 四个算法的定位差异CoSaMPCompressive Sampling Matching Pursuit每次迭代选取 ( 2K ) 个最大相关原子然后通过最小二乘估计信号并裁剪回 ( K ) 个最大分量。SPSubspace Pursuit类似但每轮选取 ( K ) 个原子并做子空间追踪收敛速度通常比 OMP 快。SAMPSparsity Adaptive Matching Pursuit不需要预先知道稀疏度 ( K )它通过逐步增加步长来逼近真实稀疏度适合实际场景中 ( K ) 未知的情况。ROMPRegularized OMP在原子选择时引入正则化分组避免选到相关性过高的原子。如果你只是做对比实验建议先跑CS_SAMP.m因为它对 ( K ) 不敏感容易出结果如果你要验证“已知稀疏度下的重构精度”那就用CS_CoSaMP.m或CS_SP.m。2.3 运行前的环境检查在 MATLAB 命令窗口执行ver确认版本R2016b 之后的版本对脚本内函数定义兼容性更好。检查当前文件夹是否已切换到cs-code.zip解压后的目录不要直接在压缩包预览里点运行。执行which CoSaMP确认路径无冲突。如果之前装过其他压缩感知工具箱比如sparsify或CVX注意ROMP这个名字可能被占用必要时用which -all ROMP排查。常见做法是新建一个run_cs_demo.m把五个文件里的测试逻辑统一到一个入口方便反复调参。3. 把 CoSaMP 和 SAMP 跑起来参数、代码与观测矩阵构造3.1 构造稀疏信号与测量矩阵的代码模板下面这段代码是我根据包内文件结构补全的通用测试框架你可以直接抄到 MATLAB 里运行。它生成一个长度 ( N256 )、稀疏度 ( K20 ) 的随机信号测量数 ( M128 )测量矩阵采用高斯随机矩阵并按列归一化。% test_cs_demo.m clear; clc; close all; N 256; % 信号长度 K 20; % 稀疏度 M 128; % 测量数 x zeros(N,1); idx randperm(N, K); x(idx) randn(K,1); % 稀疏信号非零位置随机 Phi randn(M, N); Phi Phi ./ repmat(sqrt(sum(Phi.^2,1)), M, 1); % 列归一化 y Phi * x; % 观测向量 % 调用 CoSaMP x_cosamp CoSaMP(y, Phi, K); err_cosamp norm(x - x_cosamp) / norm(x); fprintf(CoSaMP 相对误差: %.4f\n, err_cosamp); % 调用 SAMP步长设为 4 x_samp CS_SAMP(y, Phi, 4); err_samp norm(x - x_samp) / norm(x); fprintf(SAMP 相对误差: %.4f\n, err_samp);逻辑说明randperm保证非零位置不重复Phi列归一化是贪婪类算法的常见要求否则内积幅值会被列范数干扰。CoSaMP第三个参数是稀疏度 ( K )而CS_SAMP第三个参数是初始步长不是 ( K )。如果你直接给CS_SAMP传 ( K )它可能把步长设得过大导致过估计。参数怎么改N和M的比值决定压缩率一般 ( M/N ) 在 0.3 到 0.6 之间比较稳K越大重构难度越高CoSaMP 在 ( K M/4 ) 时容易崩。3.2 观测矩阵选型高斯矩阵不是唯一解包内代码默认用高斯随机矩阵但实际项目中你可能会遇到伯努利矩阵、部分傅里叶矩阵或托普利兹矩阵。高斯矩阵与大多数稀疏基不相关满足 RIP 条件的概率高但存储和计算开销大。部分傅里叶矩阵适合 MRI 和雷达成像因为物理采样过程本身就是傅里叶采样。如果你要替换测量矩阵注意保持列归一化并且检查互相关性max(max(abs(Phi*Phi - eye(M))))这个值越小越好。常见做法是先用高斯矩阵验证算法逻辑再换成特定应用的矩阵做仿真。3.3 迭代停止条件的参数调整CoSaMP 和 SP 通常以残差范数小于阈值或迭代次数达到上限作为停止条件。包内代码可能硬编码了tol 1e-6和max_iter 100。如果你发现算法提前退出但误差很大先把tol放宽到1e-3试试如果迭代次数不够改成max_iter 500。SAMP 的停止条件更复杂它比较相邻两次估计的能量差如果步长增加后能量不再显著下降就停止。你可以手动在CS_SAMP.m里找到if norm(x_new - x_old) 1e-6这类行把阈值调小会让它跑更久但精度未必提升反而可能过拟合噪声。4. 避坑与排查五个让你少熬两晚的血泪经验4.1 现象运行后报“索引超出矩阵维度”原因CoSaMP.m内部在裁剪支持集时用了x(Omega)而Omega是列向量但你的y是行向量导致Phi * y维度不匹配。解决在调用前统一y y(:);强制转成列向量。这个坑我踩过MATLAB 不会提示“方向不对”只会说索引越界。4.2 现象重构误差始终在 0.8 以上几乎没恢复原因测量矩阵没有列归一化或者稀疏信号的非零值方差太小被噪声淹没。解决检查Phi每列的 2-范数是否为 1把x(idx) randn(K,1) * 10放大信号幅度。另外确认K是否超过M/4超过就换 SAMP 或减小K。4.3 现象SAMP 迭代到一半报“步长超过信号长度”原因CS_SAMP.m里步长累加逻辑没有加边界判断当真实稀疏度接近 ( N ) 时步长会一直涨。解决在步长更新行后面加if step N/2, step N/2; end。或者直接限制最大迭代次数为ceil(N/step)。4.4 现象CoSaMP 和 CS_CoSaMP 结果不一致原因两个文件可能用了不同的最小二乘求解方式一个用pinv一个用Phi(:,Omega)\y。pinv默认截断奇异值数值稳定性更好但速度慢。解决统一用x_est zeros(N,1); x_est(Omega) Phi(:,Omega) \ y;并确保Omega去重。如果你要对比算法固定求解器比固定参数更重要。4.5 现象换用部分傅里叶矩阵后算法完全失效原因部分傅里叶矩阵的列之间相关性比高斯矩阵高贪婪类算法容易选错原子。解决改用 CoSaMP 或 SP它们对相关性容忍度稍好或者增加测量数 ( M )把压缩率降到 0.3 以下。如果必须用低压缩率考虑换 BP 类凸优化算法但那不在这个包里。5. 进阶技巧用蒙特卡洛扫参数把“玄学”变成曲线最后一章说一个我常用的验证方法不要只跑一次就下结论。压缩感知重构算法的性能对随机种子敏感单次实验的误差可能波动 20% 以上。我一般会写一个外层循环对每个 ( (K, M) ) 组合跑 100 次蒙特卡洛统计成功重构概率——成功定义为相对误差小于 ( 10^{-3} )。下面这段代码可以直接扩展test_cs_demo.m% mc_sweep.m K_list 10:5:40; M_list 80:16:160; trials 100; success_rate zeros(length(K_list), length(M_list)); for i 1:length(K_list) for j 1:length(M_list) succ 0; for t 1:trials N 256; K K_list(i); M M_list(j); x zeros(N,1); idx randperm(N, K); x(idx) randn(K,1); Phi randn(M, N); Phi Phi ./ repmat(sqrt(sum(Phi.^2,1)), M, 1); y Phi * x; x_hat CoSaMP(y, Phi, K); if norm(x - x_hat)/norm(x) 1e-3 succ succ 1; end end success_rate(i,j) succ / trials; end end imagesc(M_list, K_list, success_rate); colorbar; xlabel(测量数 M); ylabel(稀疏度 K); title(CoSaMP 成功重构概率);逻辑说明外层双重循环遍历 ( K ) 和 ( M )内层trials次独立实验。success_rate矩阵的每个元素是成功次数除以总次数。用imagesc画出热力图你能直观看到相变曲线——当 ( M ) 低于某个阈值时成功概率从 1 骤降到 0。参数怎么读如果曲线在 ( M/N 0.4 ) 附近陡降说明你的算法实现基本正确如果曲线平缓或始终低于 0.5检查测量矩阵归一化和最小二乘求解。这个技巧帮我排除了好几次“算法写错但单次结果看起来还行”的隐患。从那以后我每次拿到新的重构代码都强制先跑一遍蒙特卡洛相变图再去看具体信号的重构波形。希望帮到你。本文还有配套的精品资源点击获取

读完文章,也想定制专属网站?

尧图设计师 24 小时内与您沟通定制方案

免费获取报价 →
↑