简介:这份资源是面向信号处理、图像重建与无线通信方向的学生和研究人员整理的压缩感知重构算法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 四个算法的定位差异
CoSaMP(Compressive Sampling Matching Pursuit)每次迭代选取 ( 2K ) 个最大相关原子,然后通过最小二乘估计信号并裁剪回 ( K ) 个最大分量。SP(Subspace Pursuit)类似,但每轮选取 ( K ) 个原子并做子空间追踪,收敛速度通常比 OMP 快。SAMP(Sparsity Adaptive Matching Pursuit)不需要预先知道稀疏度 ( K ),它通过逐步增加步长来逼近真实稀疏度,适合实际场景中 ( K ) 未知的情况。ROMP(Regularized 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 里运行。它生成一个长度 ( N=256 )、稀疏度 ( K=20 ) 的随机信号,测量数 ( M=128 ),测量矩阵采用高斯随机矩阵并按列归一化。
% 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,检查测量矩阵归一化和最小二乘求解。这个技巧帮我排除了好几次“算法写错但单次结果看起来还行”的隐患。从那以后我每次拿到新的重构代码,都强制先跑一遍蒙特卡洛相变图,再去看具体信号的重构波形。希望帮到你。
本文还有配套的精品资源,点击获取