简介：这份资源面向学习压缩感知（Compressive Sensing）的学生与科研人员，提供在Matlab环境下实现CS算法的完整代码包，帮助读者理解稀疏表示、测量矩阵设计与信号重构的完整流程。压缩包共1436个文件，约10.23MB，以1230个.m函数脚本为主体，辅以65个C源码、34个mexw64与4个mexw32加速模块、13个.mat数据文件及若干txt、asc、raw样本数据，另有readme与pdf说明文档，兼顾算法实现与实验数据。核心代码fcsa涵盖稀疏基构造、随机与结构化测量矩阵、L1范数优化及ISTA/FISTA等重构算法，并可用PSNR、SSIM评估重构质量，内容预览中的地震、核磁等信号样本便于直接开展实验。目前已有139人学习下载，适合作为图像处理、无线通信与医学成像方向研究者的入门模板与调试参考。

1. 从 CScode.zip 说起：压缩感知为什么值得你花一个周末跑通

你手里有一份叫 CScode.zip 的代码包，名字里带着 compressive sensing。第一次看到这个组合，多数人的反应是：压缩感知不是十几年前就火过的理论吗，现在拿它出来还能干什么？我当初也这么想，直到在一个采样成本极高的项目里被逼着重读了一遍。压缩感知真正值钱的地方，不是“压缩”两个字，而是它允许你用远少于奈奎斯特要求的观测数，把原始信号恢复出来。换句话说，当采集一次数据的代价（时间、硬件、辐射剂量、带宽）高到肉疼时，这套东西能直接帮你省掉一大半的采样动作。CScode.zip 这类代码包，通常就是把稀疏表示、观测矩阵、重构算法这三件事打包成可跑的脚本，让你不用从公式推导开始，直接拿自己的数据试。它适合两类人：一类是信号处理、成像、通信方向的工程师，手里有真实采集数据但采样受限；另一类是做算法验证的研究者，想快速对比 OMP、CoSaMP、ISTA 这些重构算法的实际表现。这篇笔记就按“理论先立住、再动手能复现”的顺序，把压缩感知从概念讲到能跑通的代码，再讲清楚参数怎么设、坑在哪。

2. 压缩感知的三个前提：稀疏、非相干、重构算法

2.1 稀疏性到底在说什么

压缩感知能成立的第一块基石是稀疏性。一个长度为 N 的信号 x，如果在某个变换域 Ψ 下只有 K 个非零系数（K 远小于 N），我们就说它是 K-稀疏的。注意这里的关键词是“在某个变换域下”，不是原始时域。比如一张自然图像在时域上几乎处处非零，但在小波域或 DCT 域下，大部分系数接近零，只有少数大系数承载了主要能量。这就是为什么压缩感知在成像领域特别吃香——图像天然具备变换域稀疏性。

实际工程里没有绝对稀疏的信号，只有“可压缩”的信号，也就是系数按幅度排序后衰减很快。判断你的数据能不能用压缩感知，最直接的办法是把信号做一次变换，画出系数幅度排序曲线。如果曲线在前 10% 就掉到很小，那基本可以用；如果系数幅度均匀分布，那压缩感知帮不了你，别硬上。

这里有个容易翻车的地方：很多人拿原始像素值直接做稀疏性判断，发现不稀疏就放弃了。正确做法是先选变换基。常见选择是小波（图像）、DCT（块状信号）、FFT（频域稀疏信号）。CScode.zip 里一般会提供几个变换矩阵的生成函数，你需要在跑重构之前先确认自己的信号在哪个基下最稀疏。

2.2 观测矩阵为什么要满足 RIP

第二个前提是观测矩阵 Φ 和稀疏基 Ψ 之间的非相干性。压缩感知的观测过程是 y = Φx = ΦΨs，其中 y 是 M 维观测向量，M 远小于 N。要从 M 个观测里恢复 N 维信号，直觉上是不可能的，但如果 Φ 和 Ψ 不相干，信息就被“打散”了，恢复成为可能。

严格的理论保证是 RIP（有限等距性质），但 RIP 验证是 NP 难的，实际中没人真去验证。工程上的替代方案是看互相干系数 μ(Φ,Ψ)，值越小越好。随机高斯矩阵、随机伯努利矩阵和部分傅里叶矩阵是三种最常用的观测矩阵，它们和大多数固定基的互相干性都足够低。

我一般会这样选：如果信号是图像、变换基是小波，观测矩阵用随机高斯矩阵最省心，虽然硬件实现难，但仿真验证阶段够用；如果考虑硬件实现，部分傅里叶矩阵更实际，因为它对应的是随机欠采样，很多成像系统本来就能做。CScode.zip 里如果同时提供了这几种矩阵的生成脚本，建议先用高斯矩阵跑通流程，再换成你硬件能实现的矩阵做对比。

2.3 重构算法：从 OMP 到 ISTA 的选型逻辑

重构算法分三大类：贪婪算法（OMP、CoSaMP、SP）、凸优化方法（ISTA、FISTA、ADMM）、贝叶斯方法（SBL）。选哪个不是看哪个高级，而是看你的信号维度、稀疏度和实时性要求。

OMP 是最容易上手的，迭代选原子、最小二乘更新，代码短、调试直观。缺点是稀疏度 K 需要预先知道，而且 K 估计不准时效果波动大。CoSaMP 改进了这一点，每次迭代选多个原子再剪枝，对 K 的鲁棒性更好，但实现稍复杂。ISTA 属于凸优化路线，不需要知道 K，收敛有理论保证，但收敛速度慢，FISTA 加了动量加速后实用性强很多。

我的经验是：N 在几百到几千、K 在几十以内，OMP 和 CoSaMP 足够；N 上万、信号结构复杂，直接上 FISTA 或 ADMM。CScode.zip 这类包通常会把几种算法都放在一个目录下，用统一的接口调用，方便你横向对比。下面这段代码展示了一个典型的 OMP 重构流程，你可以直接拿去改：

import numpy as np def omp(y, Phi, K, tol=1e-6): """ y: 观测向量 (M,) Phi: 观测矩阵 (M, N) K: 稀疏度 tol: 残差阈值 """ M, N = Phi.shape residual = y.copy() idx_set = [] # 已选原子索引 x_hat = np.zeros(N) # 重构信号 for _ in range(K): # 计算残差与每个原子的内积，选最大 correlations = np.abs(Phi.T @ residual) new_idx = np.argmax(correlations) idx_set.append(new_idx) # 在已选原子集合上做最小二乘 Phi_sel = Phi[:, idx_set] x_sel, _, _, _ = np.linalg.lstsq(Phi_sel, y, rcond=None) x_hat[idx_set] = x_sel # 更新残差 residual = y - Phi_sel @ x_sel if np.linalg.norm(residual) < tol: break return x_hat

这段代码的逻辑很直白：每轮找和当前残差最相关的原子，加入支撑集，然后在这个支撑集上重新做最小二乘，更新残差。参数 K 控制迭代次数，设大了会过拟合噪声，设小了恢复不全。tol 是提前终止条件，残差足够小就停。实际使用时，Phi 的列最好做归一化，否则内积比较会被列范数大的原子主导，选错原子。

提示：OMP 对观测矩阵的列归一化非常敏感，跑之前先做Phi = Phi / np.linalg.norm(Phi, axis=0)，这一步能省掉你大量调试时间。

3. 把 CScode.zip 跑起来：环境、数据、参数三步走

3.1 环境准备与目录结构确认

拿到一个压缩感知代码包，第一步不是急着跑主脚本，而是先看清楚目录结构。典型的 CScode.zip 解压后会有这么几类文件：utils/放变换矩阵和观测矩阵生成函数，algorithms/放各种重构算法，demo/放示例脚本，data/放测试信号或图像。先打开 demo 脚本看一眼它依赖哪些函数、数据从哪来，比盲目python main.py高效得多。

环境方面，压缩感知代码对 Python 版本和科学计算库比较敏感。我一般用 conda 建一个干净环境，装 numpy、scipy、matplotlib 三件套就够跑大部分 demo。如果代码里用了 sklearn 的OrthogonalMatchingPursuit或者Lasso，再补一个 scikit-learn。不要一上来就装一堆深度学习框架，压缩感知的经典算法用不到。

conda create -n cs_env python=3.9 conda activate cs_env pip install numpy scipy matplotlib scikit-learn

装完之后先跑一个最小验证：生成一个稀疏信号，做一个随机高斯观测，用 OMP 重构，看能不能恢复。这一步能跑通，说明环境和代码基本兼容，再去碰真实数据。

3.2 构造你的第一个可复现实验

最小验证实验的目的是把“稀疏信号生成 → 观测 → 重构 → 评估”这条链路走通。下面这段代码可以直接复制运行，它不依赖 CScode.zip 里的任何文件，帮你先建立手感：

import numpy as np import matplotlib.pyplot as plt np.random.seed(42) N = 512 # 信号长度 K = 20 # 稀疏度 M = 128 # 观测数 # 1. 生成 K-稀疏信号 x = np.zeros(N) support = np.random.choice(N, K, replace=False) x[support] = np.random.randn(K) # 2. 生成观测矩阵（随机高斯）和观测向量 Phi = np.random.randn(M, N) / np.sqrt(M) y = Phi @ x # 3. OMP 重构（复用上一节的函数） x_hat = omp(y, Phi, K) # 4. 评估 mse = np.linalg.norm(x - x_hat) / np.linalg.norm(x) print(f"相对重构误差: {mse:.4f}") plt.figure(figsize=(10, 4)) plt.subplot(1, 2, 1) plt.stem(x, markerfmt='o', linefmt='b-', basefmt=' ') plt.title('原始信号') plt.subplot(1, 2, 2) plt.stem(x_hat, markerfmt='o', linefmt='r-', basefmt=' ') plt.title('OMP 重构') plt.tight_layout() plt.show()

这段代码里，M 和 K 的关系决定了能不能恢复。经验规则是 M ≥ 4K 左右比较稳，M 太小恢复会失败，M 太大就失去了压缩的意义。你可以把 M 从 128 往下调，观察误差什么时候开始爆炸，这个临界点就是你实际系统里观测数的最低要求。seed 固定是为了可复现，真实实验里要多跑几个 seed 取平均。

3.3 关键参数怎么设：M、K、迭代次数

参数设置是压缩感知落地最容易被忽视的环节。很多人跑完 demo 觉得效果不错，换自己的数据就崩了，问题往往出在参数没跟着数据走。

观测数 M 的选取：理论上 M = O(K log(N/K))，但常数因子取决于观测矩阵和重构算法。我的做法是从 M = 4K 开始试，逐步降低，找到误差还能接受的最小 M。如果你不知道 K，先用一个估计值，比如信号变换后保留 95% 能量的系数个数。

稀疏度 K 的选取：OMP 和 CoSaMP 需要 K 作为输入。K 设小了恢复不全，设大了会引入噪声原子。一个实用技巧是先用 ISTA 跑一遍（不需要 K），看重构系数的幅度分布，把明显大于噪声水平的系数个数作为 K 的估计。

迭代次数：OMP 的迭代次数就是 K，CoSaMP 一般设 K 到 2K，ISTA/FISTA 需要设最大迭代次数和收敛阈值。FISTA 通常几百次迭代就够，ISTA 可能要几千次。如果你发现重构结果还在缓慢变化，说明迭代不够；如果结果震荡，说明步长或正则参数有问题。

注意：不同算法的参数不能直接套用。OMP 的 K 和 ISTA 的 λ 是两套体系，换算法时参数要重新调。

4. 避坑与排查：压缩感知落地最常见的五个翻车点

4.1 重构误差大，但不知道是哪个环节的问题

现象：跑完重构，相对误差 0.5 以上，信号完全没恢复出来。

原因：可能是观测矩阵和稀疏基不相干性差，也可能是 K 估计错误，还可能是观测数 M 不够。三个环节都可能出问题，但报错信息不会告诉你。

解决：分步排查。先固定观测矩阵和 M，用真实支撑集做一次最小二乘，看能不能恢复。如果这一步误差就大，说明 M 不够或矩阵条件数太差；如果这一步好但 OMP 差，说明 K 设错了。再固定 K，换不同的观测矩阵对比，排除矩阵问题。这个分步法能帮你快速定位。

4.2 观测矩阵列没归一化导致 OMP 选错原子

现象：OMP 每次选的原子位置和真实支撑集对不上，重构结果稀疏但位置全错。

原因：OMP 用内积大小选原子，如果 Phi 的列范数差异大，范数大的列会被优先选中，即使它和残差的相关性并不高。

解决：在调用 OMP 之前对 Phi 做列归一化。这一步在仿真里容易被忽略，因为随机高斯矩阵的列范数本来就接近，但换成部分傅里叶矩阵或实测矩阵时，列范数差异会很大。

# 列归一化，放在 OMP 之前 norms = np.linalg.norm(Phi, axis=0, keepdims=True) Phi = Phi / norms # 注意：重构后要还原尺度 x_hat = x_hat / norms.flatten()

4.3 图像分块处理时块间伪影严重

现象：用压缩感知做图像重构，每个块单独处理后拼起来，块边界有明显不连续。

原因：分块独立重构忽略了块间相关性，每个块的稀疏度和噪声水平不同，重构质量不一致，拼起来就有伪影。

解决：两个方向。一是用重叠分块，重构后做加权平均，重叠区域取平均能平滑边界。二是改用整图重构，用二维小波变换代替块 DCT，让稀疏表示覆盖全局。重叠分块实现简单，适合快速验证；整图重构效果更好，但计算量大。

4.4 噪声环境下稀疏度估计偏大

现象：有噪声时，ISTA 重构出的系数比真实稀疏度多很多，很多小系数其实是噪声。

原因：ISTA 的 L1 正则会把噪声也拟合进去，λ 设小了噪声系数压不住。

解决：增大 λ，或者用去偏步骤（debiasing）：先用 ISTA 找到支撑集，再在支撑集上做最小二乘。去偏能显著改善噪声下的重构质量，代价是需要多一步计算。λ 的选取可以用交叉验证，或者观察重构误差随 λ 的变化曲线，选拐点。

4.5 实测数据上效果远差于仿真

现象：仿真里 M=4K 就能完美恢复，实测数据同样 M 却完全不行。

原因：仿真信号是理想 K-稀疏，实测信号只是可压缩，系数衰减不够快；实测噪声也不是高斯白噪声，可能有色或相关。

解决：先对实测信号做变换，看系数衰减曲线，估计有效稀疏度。如果衰减慢，压缩感知的压缩率就要降低，M 要加大。噪声方面，如果噪声有色，可以在重构前做预白化，或者用对噪声鲁棒性更好的算法（如 BPDN 代替 BP）。不要拿仿真参数直接套实测，这是最常见的翻车原因。

5. 进阶技巧：用交叉验证选正则参数，用去偏提升重构精度

5.1 交叉验证选 λ 的具体做法

ISTA/FISTA 的正则参数 λ 控制稀疏惩罚强度，选不好要么欠拟合要么过拟合。交叉验证是靠谱的办法，但压缩感知的观测数少，不能像普通回归那样留出验证集。我的做法是：把观测向量 y 随机分成训练集和验证集，用训练集重构，在验证集上算误差，选验证误差最小的 λ。

def cross_validate_lambda(y, Phi, lambdas, train_ratio=0.8): M = len(y) idx = np.random.permutation(M) n_train = int(M * train_ratio) train_idx, val_idx = idx[:n_train], idx[n_train:] best_lambda, best_err = None, np.inf for lam in lambdas: x_hat = fista(y[train_idx], Phi[train_idx, :], lam) err = np.linalg.norm(y[val_idx] - Phi[val_idx, :] @ x_hat) if err < best_err: best_err, best_lambda = err, lam return best_lambda

这段代码的核心思路是：用一部分观测做重构，用另一部分观测做验证。验证误差反映的是重构信号对未见观测的预测能力，比训练误差更能反映泛化性能。λ 的候选集一般取对数刻度，比如 1e-4 到 1e-1 之间取 10 个值。注意验证集不能太小，否则误差估计方差大，一般留 20% 左右。

5.2 去偏步骤为什么能提升精度

L1 正则的副作用是压缩大系数，导致重构信号幅度偏小。去偏的思路是：先用 L1 方法（ISTA/FISTA）找到支撑集，也就是哪些位置有非零系数，然后在这个支撑集上做无正则的最小二乘，恢复准确的幅度。

def debias(y, Phi, x_l1, threshold=1e-3): support = np.where(np.abs(x_l1) > threshold)[0] if len(support) == 0: return x_l1 Phi_sel = Phi[:, support] x_sel, _, _, _ = np.linalg.lstsq(Phi_sel, y, rcond=None) x_debiased = np.zeros_like(x_l1) x_debiased[support] = x_sel return x_debiased

threshold 控制哪些系数被认为是有意义的，设太小会把噪声系数也纳入支撑集，设太大会漏掉弱信号。我一般取 L1 重构系数最大值的 1% 到 5%。去偏之后，重构误差通常能降一个档次，尤其是对幅度动态范围大的信号。

5.3 一个我常用的验证习惯

每次跑完重构，我不会只看一个相对误差数字，而是画三张图：原始信号、重构信号、误差信号。误差信号如果看起来像噪声，说明重构基本到位；如果误差信号有结构，说明还有系统性偏差，可能是支撑集选错或参数没调好。这个习惯帮我抓过好几次隐蔽的 bug，比如观测矩阵行列对应关系搞反了，单看误差数字发现不了，看误差图就一目了然。

另外，我习惯在代码里固定随机种子，但验证时至少跑 10 个不同种子取平均。单次实验的运气成分很大，尤其是 M 接近临界值时，不同种子可能一个成功一个失败。取平均后的误差曲线才是可靠的判断依据。

希望这些经验帮到你，少走几个我当年踩过的坑。

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