简介这份资源是面向信号处理、图像重建与无线通信方向的学习者和研究人员整理的压缩感知重构算法MATLAB代码合集重点解决从少量测量值中恢复稀疏或可压缩信号时如何选择与对比不同重构策略的问题。压缩包共5个文件全部为.m脚本整体约7KB涵盖CoSaMP、SAMP、SP、ROMP等典型重构算法的实现便于在MATLAB环境中直接运行、修改参数并观察重构效果。已有149人学习关注说明其在课程实验与算法入门中具有一定参考价值。读者可借助这些脚本理解贪婪类算法的迭代流程、支撑集更新方式与残差收敛过程并横向比较不同算法在重构精度和运行效率上的差异为后续优化测量矩阵、调整稀疏度或设计新重构方案提供可复用的代码基础。1. 压缩感知重构算法选型从 cs-code.zip 里翻出能跑通的那一套手里拿到一个叫cs-code.zip的压缩感知代码包打开一看目录里躺着 OMP、CoSaMP、SAMP、IRLS 一堆.m和.py文件README 只有两行这时候最怕的不是算法难而是不知道先跑哪个、参数怎么设、跑出来为什么和论文对不上。压缩感知CS重构算法解决的核心问题就一句话用远少于奈奎斯特采样率的观测值把原始稀疏信号恢复出来。它适合做雷达成像、单像素相机、MRI 加速、宽带频谱感知、无线传感器数据回传的从业者也适合正在复现 CoSaMP、SAMP 这类经典算法、被支撑集选错和步长卡住的新手。这篇笔记不假设你手里有某份特定源码而是按这类代码包最常见的组织方式把选型、跑通、调参、排错整条链路讲清楚让你拿到任何一个 CS 重构代码包都能快速判断值不值得投入、怎么改到自己的场景里。2. 压缩感知重构的数学骨架与算法选型为什么 OMP 快而 CoSaMP 稳2.1 从 yAx 到稀疏恢复三个矩阵维度先对齐压缩感知的观测模型写成y A x n其中x是长度 N 的原始信号A是 M×N 的观测矩阵也叫测量矩阵、感知矩阵y是长度 M 的观测向量n是噪声。稀疏性假设是x里非零元素个数 K 远小于 N通常 K M N。重构算法要做的就是从已知的A和y里反推出x这是一个欠定问题没有稀疏约束时解不唯一加上 L0 或 L1 约束才有唯一解。落到代码里第一件要确认的事是矩阵维度方向。很多代码包默认A是 M×N、x是 N×1、y是 M×1但有些包为了 MATLAB 的列优先习惯把信号写成行向量转置没对齐就会报维度不匹配或者更坑的——不报错但结果全错。我一般会在跑主脚本前先插一段维度检查import numpy as np # 假设从代码包里加载观测矩阵和观测向量 A np.load(A_matrix.npy) # 期望形状 (M, N) y np.load(y_vector.npy) # 期望形状 (M,) 或 (M, 1) x_true np.load(x_true.npy) # 期望形状 (N,) 或 (N, 1) M, N A.shape print(fA: {A.shape}, y: {y.shape}, x_true: {x_true.shape}) # 统一压成列向量避免后续广播出错 y y.reshape(-1, 1) x_true x_true.reshape(-1, 1) assert y.shape[0] M, fy 长度 {y.shape[0]} 与 A 行数 {M} 不一致 assert x_true.shape[0] N, fx_true 长度 {x_true.shape[0]} 与 A 列数 {N} 不一致 # 稀疏度非零元素个数阈值按信号量级取 K np.sum(np.abs(x_true) 1e-6) print(f稀疏度 K {K}, 压缩比 M/N {M/N:.3f})这段代码做三件事打印三个对象的形状、把向量统一成列向量、用 1e-6 的阈值估计稀疏度 K。K 这个值非常关键OMP、CoSaMP、SAMP 里凡是带「已知稀疏度」的算法都要用它估错了后面全崩。阈值不要死记 1e-6要看信号本身的量级如果x_true元素普遍在 1e-3 量级阈值就该放到 1e-8 甚至更小否则会把真实非零值当成零。2.2 贪心类、阈值类、凸优化类三类算法的适用边界CS 重构算法大致分三派选型时先看你的场景属于哪一类。贪心类包括 OMP、CoSaMP、SAMP、Subspace Pursuit。它们的思路是每次迭代挑出和残差最相关的若干列加入支撑集再最小二乘更新。OMP 每次只加一列速度快但需要已知 KCoSaMP 每次加 2K 列再剪回 K 列抗噪更好但单次迭代开销大SAMP 不需要已知 K靠步长逐步扩大支撑集适合 K 未知的场景。贪心类的共同短板是当观测矩阵列相关性高时容易选错原子一旦选错很难回头。阈值类包括 IHT、Hard Thresholding、Soft Thresholding。每步做梯度下降再阈值化实现简单适合大规模问题但收敛对步长敏感步长设大了直接发散。凸优化类把 L0 松弛成 L1用 ISTA、FISTA、ADMM 求解。重构精度通常最好代价是迭代次数多、计算量大适合对精度要求高、可以离线跑的场合。选型建议按这个顺序问自己K 已知吗已知就 OMP 或 CoSaMP 起步未知就 SAMP 或 IHT。观测矩阵列相关性高吗高就优先凸优化类。实时性要求高吗高就贪心类低就凸优化类。下面这张表是我在多个项目里总结的对照参数范围是常见取值不是绝对标准。算法需要已知 K典型迭代次数抗噪表现适用场景OMP是K 到 2K一般快速原型、K 已知CoSaMP是5 到 10较好中等规模、有噪声SAMP否10 到 30较好K 未知、稀疏度变化IHT是50 到 200一般大规模、简单实现FISTA否100 到 500好高精度、离线2.3 观测矩阵怎么选随机高斯为什么是默认答案观测矩阵A决定了重构的下界。常见选择有随机高斯矩阵、随机伯努利矩阵、部分傅里叶矩阵、托普利兹矩阵。随机高斯矩阵的优点是和绝大多数稀疏基不相关满足 RIP 条件的概率高代码包里默认用它最多。部分傅里叶矩阵适合有快速变换的场景比如 MRI 里 K 空间采样天然就是傅里叶观测计算快但列相关性比高斯高。如果你要自己生成观测矩阵常见做法是import numpy as np def generate_gaussian_matrix(M, N, seed42): 生成 MxN 随机高斯观测矩阵列归一化 rng np.random.default_rng(seed) A rng.standard_normal((M, N)) # 列归一化每列 L2 范数归一到 1避免量级差异影响内积比较 norms np.linalg.norm(A, axis0, keepdimsTrue) A A / norms return A M, N 256, 1024 A generate_gaussian_matrix(M, N) print(f观测矩阵形状: {A.shape}, 列范数示例: {np.linalg.norm(A[:, 0]):.4f})列归一化这一步很多人会漏。贪心算法靠内积A^T r挑原子如果各列范数不一致范数大的列天然占便宜选出来的支撑集就偏了。归一化之后每列公平竞争重构成功率明显提升。种子seed固定是为了结果可复现调试阶段一定要固定上线再考虑换随机源。3. 把 cs-code.zip 跑起来从目录结构到第一个重构结果3.1 代码包目录的典型结构与入口定位这类 CS 重构代码包不管作者怎么组织通常逃不出这几种文件算法实现文件omp.m、cosamp.py、samp.py、工具函数generate_A.m、compute_snr.m、演示脚本demo_omp.m、test_all.py、数据文件.mat、.npy。拿到包先别急着跑用命令行把结构列出来# 列出目录树只看两层过滤掉版本控制目录 find . -maxdepth 2 -type f \( -name *.m -o -name *.py -o -name *.mat -o -name *.npy \) | sort # 统计各类文件数量快速判断包的重心 find . -name *.m | wc -l find . -name *.py | wc -l如果.m文件占绝大多数说明是 MATLAB 包需要 MATLAB 或 Octave 环境如果.py多就是 Python 包检查有没有requirements.txt或setup.py。入口脚本通常名字里带demo、main、test、run优先找这些。找不到入口就找 README 或注释最密集的那个文件作者一般会在里面写调用示例。3.2 最小可复现实验合成稀疏信号 高斯矩阵 OMP不管代码包自带什么数据我都会先跑一个自己构造的最小实验把算法本身和代码包的数据加载逻辑解耦。这样如果自带数据跑不通能立刻判断是算法问题还是数据问题。import numpy as np def omp(A, y, K, tol1e-6, max_iterNone): 正交匹配追踪 A: 观测矩阵 (M, N) y: 观测向量 (M, 1) K: 稀疏度 tol: 残差范数停止阈值 M, N A.shape if max_iter is None: max_iter K * 3 # 防止 K 估错时死循环 residual y.copy() support [] # 已选原子索引 x_hat np.zeros((N, 1)) for it in range(max_iter): # 计算残差与所有原子的内积取绝对值最大者 correlations np.abs(A.T residual).flatten() # 排除已选索引避免重复选同一列 correlations[support] -1 new_idx int(np.argmax(correlations)) support.append(new_idx) # 在支撑集上做最小二乘 A_sub A[:, support] x_sub, _, _, _ np.linalg.lstsq(A_sub, y, rcondNone) residual y - A_sub x_sub # 残差足够小就提前退出 if np.linalg.norm(residual) tol: break x_hat[support] x_sub return x_hat, support, it 1 # 构造实验 np.random.seed(0) N, M, K 512, 128, 20 A np.random.randn(M, N) A A / np.linalg.norm(A, axis0, keepdimsTrue) x_true np.zeros((N, 1)) idx np.random.choice(N, K, replaceFalse) x_true[idx] np.random.randn(K, 1) y A x_true x_hat, support, n_iter omp(A, y, K) # 重构误差相对 L2 误差 err np.linalg.norm(x_hat - x_true) / np.linalg.norm(x_true) print(f迭代次数: {n_iter}, 相对误差: {err:.4e}, 支撑集命中: {len(set(support) set(idx))}/{K})这段 OMP 实现里有几个关键点。correlations[support] -1是防止重复选列少了这句算法会在同一列上反复打转。max_iter K * 3是保险丝K 估小了还能多迭代几次补救K 估大了也不至于跑到天荒地老。np.linalg.lstsq的rcondNone在新版 NumPy 里必须显式写否则报警告。最后那个「支撑集命中」指标比误差更直观命中数低于 K 的 80% 就说明选原子出了问题得回头查观测矩阵或 K 的估计。3.3 从 OMP 切到 CoSaMP 和 SAMP接口差异与参数映射同一个代码包里换算法最容易翻车的地方是参数名和接口不一致。OMP 要 KCoSaMP 也要 KSAMP 不要 K 但要步长step_size和最大支撑集上限。下面是把上面实验切到 CoSaMP 和 SAMP 的调用示意参数含义逐个说清。def cosamp(A, y, K, max_iter50, tol1e-6): 压缩采样匹配追踪 K: 稀疏度每轮保留 2K 个候选再剪回 K 个 M, N A.shape residual y.copy() support np.array([], dtypeint) x_hat np.zeros((N, 1)) for it in range(max_iter): # 步骤1找 2K 个最相关原子 corr np.abs(A.T residual).flatten() cand np.argsort(corr)[::-1][:2*K] # 步骤2合并当前支撑集与候选 union np.union1d(support, cand) # 步骤3在并集上最小二乘 A_sub A[:, union] x_sub, _, _, _ np.linalg.lstsq(A_sub, y, rcondNone) # 步骤4剪枝保留幅值最大的 K 个 keep np.argsort(np.abs(x_sub.flatten()))[::-1][:K] support union[keep] # 步骤5更新残差 A_sel A[:, support] x_sel, _, _, _ np.linalg.lstsq(A_sel, y, rcondNone) residual y - A_sel x_sel if np.linalg.norm(residual) tol: break x_hat[support] x_sel return x_hat, support, it 1 def samp(A, y, step_size4, max_supportNone, max_iter100, tol1e-6): 稀疏度自适应匹配追踪 step_size: 每轮支撑集扩大的列数 max_support: 支撑集上限默认 M/2 M, N A.shape if max_support is None: max_support M // 2 residual y.copy() support np.array([], dtypeint) stage 0 x_hat np.zeros((N, 1)) while len(support) max_support and stage max_iter: # 当前阶段的目标支撑集大小 target_size min((stage 1) * step_size, max_support) corr np.abs(A.T residual).flatten() cand np.argsort(corr)[::-1][:target_size] union np.union1d(support, cand) A_sub A[:, union] x_sub, _, _, _ np.linalg.lstsq(A_sub, y, rcondNone) keep np.argsort(np.abs(x_sub.flatten()))[::-1][:target_size] new_support union[keep] A_sel A[:, new_support] x_sel, _, _, _ np.linalg.lstsq(A_sel, y, rcondNone) new_residual y - A_sel x_sel # 残差没明显下降就进入下一阶段扩大步长 if np.linalg.norm(new_residual) np.linalg.norm(residual) * 0.99: stage 1 else: support new_support residual new_residual x_hat np.zeros((N, 1)) x_hat[support] x_sel if np.linalg.norm(residual) tol: break return x_hat, support, stage 1CoSaMP 和 OMP 的核心差别在「先扩后剪」每轮先选 2K 个候选和已有支撑集合并最小二乘之后只留 K 个。这个机制让 CoSaMP 有回头路选错的原子下一轮可能被剪掉所以抗噪比 OMP 好。max_iter50是经验值K 在几十的量级时通常 10 轮内收敛。SAMP 的关键参数是step_size。设小了阶段推进慢、迭代次数多设大了支撑集一次扩太多容易把不相关的原子也拉进来。常见做法是从 2 到 8 之间试稀疏度大就取大一点。max_support M // 2是理论上限支撑集超过 M/2 时最小二乘已经过拟合重构没有意义。残差下降判断用0.99这个系数意思是残差没降到上一轮的 99% 以下就认为当前阶段榨干了进入下一阶段。这个系数太松会导致阶段推进过快太紧会卡在原地0.95 到 0.99 之间调。4. 参数调优与性能验证K 估错、步长设偏、噪声上来之后怎么办4.1 稀疏度 K 估错的三种表现与补救K 是贪心类算法的命门。K 估小了的典型表现是重构误差降到一个平台就下不去支撑集命中数明显小于真实 K残差范数停在某个非零值。补救办法是跑一遍 SAMP 或者把 OMP 的max_iter放大到K_est * 5看误差能不能继续降。K 估大了的表现是支撑集里混入不相关原子重构信号出现很多小幅度毛刺相对误差反而比 K 准确时更大。这时候要做的是在 OMP 每轮迭代后记录残差画残差随迭代次数的曲线找到残差下降明显变缓的拐点那个位置的迭代次数就是更合理的 K。# 残差曲线诊断跑 OMP 并记录每轮残差 def omp_with_history(A, y, K): M, N A.shape residual y.copy() support [] history [] for it in range(K * 5): corr np.abs(A.T residual).flatten() corr[support] -1 idx int(np.argmax(corr)) support.append(idx) A_sub A[:, support] x_sub, _, _, _ np.linalg.lstsq(A_sub, y, rcondNone) residual y - A_sub x_sub history.append(np.linalg.norm(residual)) return history hist omp_with_history(A, y, K) # 找残差下降幅度小于前一轮 10% 的位置 for i in range(1, len(hist)): drop (hist[i-1] - hist[i]) / (hist[i-1] 1e-12) if drop 0.1: print(f建议 K 不超过 {i1}此处残差下降仅 {drop:.2%}) break这段诊断代码不依赖真实 K靠残差下降速率反推。drop 0.1是经验阈值意思是这一轮加进来的原子对残差几乎没有贡献再加下去就是过拟合。4.2 步长与迭代上限SAMP 和 IHT 的敏感参数SAMP 的step_size和 IHT 的步长mu是两类最容易设偏的参数。SAMP 步长偏小的表现是迭代次数暴涨、每阶段只扩几列、跑到max_iter还没收敛偏大的表现是支撑集一次扩过头、重构误差在某个阶段突然跳高。我一般先用step_size4跑一遍看收敛时的阶段数和最终支撑集大小如果阶段数超过 15 就说明步长偏小调到 8 再试。IHT 的步长mu必须满足mu 1 / ||A||_2^2其中||A||_2是观测矩阵的谱范数。设大了直接发散残差越来越大设小了收敛慢。计算谱范数用np.linalg.norm(A, 2)然后取mu 0.9 / (norm_A ** 2)作为起点。def iht(A, y, K, muNone, max_iter200, tol1e-6): 迭代硬阈值 mu: 步长默认按谱范数自动设置 M, N A.shape if mu is None: norm_A np.linalg.norm(A, 2) mu 0.9 / (norm_A ** 2) x_hat np.zeros((N, 1)) for it in range(max_iter): # 梯度下降步 grad A.T (A x_hat - y) x_new x_hat - mu * grad # 硬阈值只保留幅值最大的 K 个 idx np.argsort(np.abs(x_new.flatten()))[::-1][:K] x_hat np.zeros((N, 1)) x_hat[idx] x_new[idx] if np.linalg.norm(A x_hat - y) tol: break return x_hat, it 1mu 0.9 / (norm_A ** 2)里的 0.9 是安全系数留 10% 余量防止边界发散。如果跑出来残差震荡不降先把mu乘 0.5 再试。4.3 噪声水平与重构成功率SNR 扫描实验怎么做验证一个 CS 重构算法值不值得用最直接的办法是扫 SNR看重构成功率随 SNR 的变化曲线。成功率定义成相对误差小于某个阈值比如 1e-2的实验次数占比。下面这段代码做 SNR 从 10dB 到 40dB 的扫描每个 SNR 跑 50 次蒙特卡洛。def run_snr_sweep(A, N, K, snr_list, n_trials50, algoomp): results {} for snr_db in snr_list: success 0 for t in range(n_trials): rng np.random.default_rng(t) x_true np.zeros((N, 1)) idx rng.choice(N, K, replaceFalse) x_true[idx] rng.standard_normal((K, 1)) y_clean A x_true # 按 SNR 加高斯噪声 signal_power np.mean(y_clean ** 2) noise_power signal_power / (10 ** (snr_db / 10)) noise rng.standard_normal(y_clean.shape) * np.sqrt(noise_power) y y_clean noise if algo omp: x_hat, _, _ omp(A, y, K) elif algo cosamp: x_hat, _, _ cosamp(A, y, K) err np.linalg.norm(x_hat - x_true) / np.linalg.norm(x_true) if err 1e-2: success 1 results[snr_db] success / n_trials return results snr_list [10, 15, 20, 25, 30, 35, 40] res_omp run_snr_sweep(A, N, K, snr_list, algoomp) res_cosamp run_snr_sweep(A, N, K, snr_list, algocosamp) for snr in snr_list: print(fSNR{snr}dB OMP{res_omp[snr]:.2f} CoSaMP{res_cosamp[snr]:.2f})n_trials50是精度和耗时的折中要更平滑的曲线就加到 200。err 1e-2这个成功率阈值按你的应用定成像类可能要求 1e-3通信类 1e-1 就够。跑完这条曲线你就能回答「这个算法在我的噪声条件下能不能用」这个问题而不是凭感觉。5. 避坑与排查CS 重构代码包里最容易翻车的 5 个地方5.1 现象重构误差始终在 0.5 以上怎么调 K 都没用原因观测矩阵没有列归一化或者归一化方式不一致。有些代码包在生成A时归一化了但加载保存过的A时又做了一次导致列范数偏离 1 很远。贪心算法靠内积选原子列范数不齐直接导致选错。解决加载A后立刻检查列范数分布np.linalg.norm(A, axis0)的均值和标准差均值应该在 1 附近标准差小于 0.01。不满足就重新归一化并且确认后续所有用到A的地方都用同一份归一化后的矩阵。5.2 现象SAMP 跑出来支撑集大小远超真实 K误差比 OMP 还大原因step_size设得太大或者残差下降判断的系数太松导致阶段推进过快支撑集一次扩过头。另一个常见原因是max_support设成了 N 而不是 M/2支撑集可以无限扩。解决先把step_size降到 2max_support显式设成M // 2残差下降系数从 0.99 收紧到 0.95。跑一遍看收敛时的支撑集大小如果还是偏大检查观测矩阵的列相关性相关性高时 SAMP 本身就不适合换 CoSaMP 或凸优化类。5.3 现象MATLAB 代码包在 Octave 里跑报错提示函数未定义原因代码包用了 MATLAB 工具箱函数比如lsqr、spgl1、cvxOctave 没有或者接口不同。另一个原因是用了arguments块做参数校验Octave 版本低不支持。解决先看报错函数名如果是工具箱函数找 Octave 的替代实现或者换成基础运算符。lsqr可以用A \ y或pinv(A) * y替代精度略降但能跑。如果是语法不兼容把arguments块改成手动nargin判断。实在改不动就装 MATLAB 或找纯 Python 版本。5.4 现象小规模实验成功放大到 N4096 后内存爆掉原因观测矩阵A是 M×N 稠密矩阵N4096、M1024 时单精度占 16MB双精度 32MB本身不大但最小二乘里A_sub.T A_sub或者np.linalg.lstsq内部会复制多份加上蒙特卡洛循环里反复分配内存峰值可能是矩阵本身的几十倍。解决把A转成float32最小二乘用scipy.linalg.lstsq并指定lapack_drivergelsy减少内存复制。蒙特卡洛循环里复用预分配的数组不要每轮np.zeros。如果还爆把A存成memmap或者用矩阵-free 的算子形式每次只算A x和A.T y不显式存A。5.5 现象同一份代码两次运行结果不一样误差差很多原因随机种子没固定。观测矩阵生成、稀疏信号构造、噪声添加任何一处用了全局随机状态且没设种子结果就不可复现。另一个隐蔽原因是多线程 BLAS 的归约顺序不确定导致最小二乘结果有微小差异在迭代算法里被放大。解决所有随机源用独立的np.random.default_rng(seed)不要用全局np.random.seed。调试阶段设OMP_NUM_THREADS1关掉 BLAS 多线程确认结果稳定后再开。如果关了多线程还不稳定检查有没有用到set或dict的遍历顺序Python 3.7 之后字典有序但集合仍然无序支撑集用list或np.array存。6. 进阶技巧用支撑集重合度和残差曲线判断算法是否真的收敛跑通重构只是第一步判断算法是不是真的收敛、结果可不可信靠单一误差指标不够。我习惯同时看两个东西支撑集重合度和残差下降曲线。支撑集重合度是重构支撑集和真实支撑集的 Jaccard 系数|S_hat ∩ S_true| / |S_hat ∪ S_true|。这个指标比误差更敏感误差 0.01 时重合度可能只有 0.7说明重构信号数值上接近但支撑集选偏了换一组观测可能就崩。重合度低于 0.9 就要警惕低于 0.7 基本不可用。def support_jaccard(support_hat, support_true): s_hat set(support_hat) s_true set(support_true) inter len(s_hat s_true) union len(s_hat | s_true) return inter / union if union 0 else 0.0 # 在蒙特卡洛里同时统计误差和重合度 jaccard_scores [] for t in range(50): rng np.random.default_rng(t) x_true np.zeros((N, 1)) idx rng.choice(N, K, replaceFalse) x_true[idx] rng.standard_normal((K, 1)) y A x_true x_hat, support_hat, _ omp(A, y, K) jaccard_scores.append(support_jaccard(support_hat, idx)) print(f平均 Jaccard: {np.mean(jaccard_scores):.3f}, 最低: {np.min(jaccard_scores):.3f})残差曲线看的是下降形态。健康的收敛曲线是指数下降每轮残差乘一个小于 1 的因子。如果曲线出现平台期再突然下降说明算法在某个阶段卡住了可能是支撑集选错后靠后续迭代补救。如果曲线震荡说明步长或 K 设得不对。把残差取对数画出来斜率应该基本恒定斜率突然变缓的位置就是该停的地方。这两个指标配合使用能帮你在没有真实信号可比的场景下判断重构结果可不可信。我现在的习惯是任何 CS 重构实验误差、重合度、残差曲线三个一起看缺一个都不放心。这套流程在雷达和频谱感知项目里救过我好几次避免了把选错支撑集的结果当成有效重构交出去。希望帮到你。本文还有配套的精品资源点击获取