简介：压缩感知（Compressed Sensing, CS）的 MATLAB 代码包，面向信号处理、数据采集、医学成像等领域的工程师、研究者与学生，聚焦多个正弦信号的随机欠采样与重构问题。包内实现正交匹配追踪（OMP）和 SPGL1 两种主流恢复算法：OMP 通过迭代挑选基原子逼近稀疏解，计算直观，适合快速验证；SPGL1 基于 L1 范数最小化，在噪声干扰或欠采样率较低时重构更稳定，二者互补，便于对比学习，也适合课程实验与算法验证。压缩包共 36 个文件，以 25 个 M 文件为核心，另有 C/H 源码和 mex 编译文件辅助运行，并附 README、ChangeLog 等说明文档，整体约 67KB，轻量易部署。已有 1717 人学习下载。借助示例脚本，可完整跑通随机欠采样、稀疏重构与误差评估流程，理解 CS 的采样和恢复机理，也可基于这些代码替换信号模型、调整参数，适配通信、雷达或成像等实际工程场景。

1. 压缩感知恢复多个正弦信号：随机欠采样为什么能行，Matlab里怎么做

很多人第一次听说压缩感知（Compressed Sensing, CS），是冲着“用远低于奈奎斯特率的采样也能恢复信号”这句断言来的。手里正好有一批振动、心电或射频信号，想赶紧在Matlab里验证一把，结果打开编辑器又不知道从哪里下手。本文针对这个最经典的起点场景——多个正弦信号的线性叠加，讲清楚随机欠采样为什么有效、OMP恢复算法怎么写、参数怎么调，以及恢复失败时先怀疑谁。目标很直接：让你拿到一段能直接跑的Matlab代码，看懂每一行在干什么，并且能自己改频率、改采样点数做实验。适合刚接触CS、准备课程设计，或者想评估压缩感知在自己信号上是否可行的工程师。

2. 稀疏表示与观测矩阵：正弦信号用少采样点恢复的前提条件是什么

压缩感知不是万能采样器。它能把采样点数压下来，前提是信号在某个变换域里有稀疏表示，而且观测方式要和这个变换域“错开”。搞懂这两个前提，再去调代码才不会一头雾水。

2.1 正弦信号在频域只有几条谱线：稀疏表示是CS的前提

一个由三个正弦叠加的实信号x(t)，如果走一遍N点FFT，结果里真正显著的谱线只有六条：正频率三条、负频率三条，其余位置接近零。换句话说，信号在时域看起来铺满整个时间轴，但在傅里叶基下只占据了少数字典元素。这就是压缩感知模型的典型场景：x = Ψs，Ψ取傅里叶基，s中非零项极少。

这个“极少”是相对N而言的。N=512时，稀疏度K只有3或6，非零占比不到2%。压缩感知理论给出的恢复边界，几乎全部由K决定，而不是由信号带宽直接决定。所以测试CS算法时，第一个动作永远是确认目标信号的稀疏度，也就是“在哪个域里稀疏、大概有几个非零系数”。正弦信号是CS天然的标准考题，因为它的稀疏性一目了然，恢复失败时便于排查是稀疏假设出问题，还是算法实现出问题。

有一点容易误解：稀疏度不等于正弦分量的个数。如果信号是三个实数余弦，频域有六个非零位置；如果用了复数指数形式的字典，则可以只数正频率的三条。本文后续代码使用傅里叶字典恢复频谱，s_hat里的非零位置会对应正负频率成对出现，画图时看正频率一侧即可。这个细节不影响恢复算法的运行，但会影响你数稀疏度。

2.2 随机欠采样怎么把欠定方程变成能解：观测矩阵与RIP的直观理解

均匀网格上有N个采样位置，现在只取其中M个，数学上就是用观测矩阵Φ抽出行，得到一个欠定方程组y = Φx。M < N，未知数个数大于方程个数，按线性代数的常识根本解不出来。CS的关键突破在于：如果x是稀疏的，并且Φ在“稀疏向量集合”上满足受限等距性质（RIP），那么欠定方程反而存在唯一解，而且这个解就是原始信号。

RIP的通俗理解是：感知矩阵ΦΨ（观测矩阵和稀疏字典合并后的整体）在任意K稀疏向量上近似保持向量长度。换句话说，稀疏向量经过这个矩阵变换后不会被“压扁”或“拉长”太多。这样，两个不同的K稀疏信号，经过观测后得到的y也不同，不会发生碰撞。严格验证RIP是NP难题，但随机矩阵有很好的理论保证：高斯随机矩阵、部分傅里叶矩阵都在高概率下满足RIP。这正是“随机欠采样”这个名字的来源。

很多初学者会问：均匀间隔抽行不也是欠采样吗？行不行？答案是几乎肯定不行。均匀抽取会让奈奎斯特采样定理意义上的频谱混叠集中爆发——原始谱线镜像折叠到低频处，和真实谱线叠在一起，稀疏性直接被破坏。随机抽取则不同，混叠能量像噪声一样平摊到整个频谱，真实谱线仍然作为大峰值冒出来，贪心算法就有机会逐步找到它们。这是随机欠采样相比均匀欠采样的本质优势。

2.3 恢复算法的选型：OMP、基追踪与阈值迭代的取舍

稀疏恢复算法大体可分为贪心类、凸优化类和迭代阈值类。OMP是贪心类代表，每次选一个与残差相关性最强的原子，然后用最小二乘更新系数，迭代K步。它的实现成本低，跑得快，在稀疏度小、信噪比高的场景下表现可靠，这也是本文选择OMP作为主算法的原因。

基追踪（Basis Pursuit）把l0问题松弛成l1凸优化，对字典原子高度相关的情况更耐受，但需要CVX或l1-magic这类外部优化工具，部署重一点。迭代硬阈值（IHT）的迭代格式简单，适合超大规模问题，但收敛慢，参数敏感性较高。三者的取舍可以简化为：先跑OMP，如果发现两个频率靠太近恢复不稳定，再上基追踪；信号维度极大时考虑IHT。

3. Matlab代码实现：生成多正弦、随机欠采样、OMP恢复的一次完整跑通

这一章给的是完整可运行的最小方案。先把信号造出来，再随机取样本，最后用OMP从欠采样数据里恢复频谱。如果你只是想验证CS能不能用在正弦信号上，把这一章跑通就够用了。

3.1 生成多正弦测试信号：频率、幅度与相位怎么设

测试信号要能暴露算法问题，不能太理想。常见的反面做法是取两个整数频率，比如50Hz和120Hz，恰好落在DFT频点上，恢复出来完美无缺，一换到真实数据立刻翻车。更合理的做法是把频率设成非整数，比如50.3Hz，让真实频率偏离DFT网格，这时候算法必须处理字典失配问题。

fs = 1024; % 等效采样率，奈奎斯特网格基准 N = 512; % 均匀网格点数，时长约 0.5 秒 t = (0:N-1).' / fs; % 时间列向量 % 三个正弦分量的参数：频率、幅度、初相位 freqs = [50.3, 120.7, 240.2]; % 故意取非整数频率，测试网格失配 amps = [1.0, 0.5, 0.3]; % 幅度差异拉开，便于观察恢复效果 phases = [0.1, 0.8, 2.0]; % 合成信号，N×1 列向量 x = zeros(N, 1); for k = 1:numel(freqs) x = x + amps(k) * cos(2*pi*freqs(k)*t + phases(k)); end

这段代码的核心是把信号生成和后续观测解耦。t用列向量，freqs用行向量，外积时可以让fortran风格的内存布局更清晰，也避免Matlab隐式扩展带来的维度混乱。相位的取值没有特殊含义，只是为了不让三个正弦在某个初始时刻恰好同相或反相，造成信号波形的偶然特性影响实验结果。

如果后续要加噪声，比较自然的写法是x = x + sigma * randn(N,1)。sigma取值可以和最大幅度关联，比如sigma = 0.05 * max(abs(x))，这样信噪比大约在26dB上下，能清楚看出OMP在噪声下的行为。

3.2 随机欠采样的Matlab实现：从均匀网格里挑点的两种写法

随机欠采样的本质，是在N个均匀网格点里无放回地抽取M个索引。Matlab里一句randperm(N, M)就足够，但有两个细节值得注意：一是要固定随机种子，否则每次跑出来的观测矩阵都不同，实验结果无法对比；二是索引要排序，后续构造观测矩阵时行的顺序是乱的，不排序会给自己找麻烦。

rng(42); % 固定种子，保证实验可复现 M = 128; % 观测样本数，约为 N 的 1/4 idx = sort(randperm(N, M).'); % 从 1:N 中无放回抽 M 个并排序 y = x(idx); % 实际观测值，M×1 列向量 % 查看采样位置分布，确认没有明显成段空洞 figure; plot(t, x, 'b-', t(idx), y, 'ro'); xlabel('时间/s'); ylabel('幅度'); legend('原始信号', '随机采样点');

另一种常见写法是用伯努利掩码：mask = rand(N,1) < p，其中p是采样比例。这种写法更接近传感矩阵理论中的“随机行选取”模型，但得到的采样点数M本身是随机的，不便于控制实验变量。实际调试时推荐用randperm方案，因为你确切知道M是多少，能准确计算采样率和恢复边界的关系。

3.3 OMP恢复核心函数：原子选择与残差更新

OMP的每一步都在做三件事：找和当前残差最匹配的原子，把它的索引加入支撑集，用最小二乘法重新拟合所有已选原子并计算新残差。重复K次，残差就逐渐被拆解到选中的原子方向上。

function s_hat = omp_recover(A, y, K, tol) % A : M×N 字典矩阵，每一列对应一个频域原子 % y : M×1 观测向量 % K : 最大迭代次数，即稀疏度上限 % tol : 残差阈值，残差范数低于该值则提前终止 [M, N] = size(A); % 原子归一化，保证相关比较时不受能量差异影响 A_norm = sqrt(sum(abs(A).^2, 1)); A_unit = A ./ A_norm; res = y; support = []; % 支撑集，记录已选原子索引 for iter = 1:K % 1. 计算每个原子与残差的内积，找最匹配原子 corr = A_unit' * res; % 复共轭转置，已包含共轭运算 [~, pos] = max(abs(corr)); % 防止同一个原子被重复选中 if ismember(pos, support) break; end support(end+1) = pos; % 2. 用最小二乘拟合当前支撑集上的系数 coef = A(:, support) \ y; % 3. 更新残差 res = y - A(:, support) * coef; % 残差足够小就提前停下来，避免继续选噪声 if norm(res) < tol break; end end % 按最终支撑集重新求解一次，得到稀疏频谱 s_hat = zeros(N, 1); s_hat(support) = A(:, support) \ y; end

这里两个容易写错的地方：一是相关系数计算必须用A_unit'而不是A_unit.'，前者是共轭转置，对复字典来说这才是内积；二是支撑集上的最小二乘必须用原始A，不能用归一化后的A_unit，否则恢复出的频谱幅度是按单位原子标定的，不是信号的真实幅度。tol参数在无噪声时可以设到1e-6，有噪声时应该和噪声水平匹配，否则OMP会把噪声也当信号选进来。

3.4 主脚本串起来：跑通一次全流程

字典矩阵在这一步构造。信号是N点均匀网格，频率网格也取N点DFT的频点：0到fs/2正频率各点，再加上负频率对应点。观测矩阵就是在这些频率原子上，只保留采样时刻idx对应的行。

% 频率网格：N 点 DFT 的频率刻度，单位 Hz f_grid = (0:N-1) * fs / N; % 部分傅里叶字典：行 = 采样时刻，列 = 频率原子 t_obs = t(idx); A = exp(1i * 2 * pi * t_obs * f_grid) / sqrt(N); % 调用 OMP 恢复频谱 K = 8; % 稀疏度上限，真实分量 6 条谱线，留余量 tol = 1e-6; s_hat = omp_recover(A, y, K, tol); % 打印最大的三个频谱峰对应频率 [mag, pos] = sort(abs(s_hat), 'descend'); fprintf('恢复的主要频率:\n'); for kk = 1:min(3, numel(mag)) fprintf(' f = %.3f Hz, 幅度 = %.3f\n', f_grid(pos(kk)), mag(kk)); end % 重构时域信号并计算误差 x_recon = real(A * s_hat); recon_err = norm(x_recon - x) / norm(x); fprintf('重构相对误差 = %.3e\n', recon_err);

跑完这段，大概率打印出接近50Hz、120Hz、240Hz的三个频率，重构误差在1e-6量级。但细心的人会发现打印的频率可能是50.000、120.000、240.000，而不是真实的50.3、120.7、240.2。原因是OMP只能在给定的频率网格上选原子，真实频率落在网格之间时，它只能选最近的网格点。这个现象就是第四章要讲的字典失配问题。

4. 三个必调参数：稀疏度、采样率、频率网格怎么设才能恢复得准

OMP跑通只是第一步。换成自己的信号时，三个参数几乎是必然要动的：稀疏度K、观测样本数M、频率网格分辨率。每个参数设错都会带来特定症状，理解症状才能对症下药。

4.1 稀疏度K对恢复的影响：设小则遗漏频率，设大则引入伪峰

K是OMP的最大迭代次数，也就是算法最多选多少个原子。K设小了，真实谱线还没选完，残差就被扔掉了，结果会漏掉某个频率分量；K设大了，算法扫完了真实分量之后不一定会停下来，而是继续在残差里找相关性最高的原子，把噪声或网格泄漏当成信号选进来，频谱上出现伪峰。

想确认K是否够用，可以不看恢复出的幅度，改看残差的变化轨迹。每次OMP迭代后残差范数应该断崖式下降，到真实分量选完时出现明显拐点，之后再迭代残差下降非常缓慢。在K未知时，扫描K=1到K=20，画残差范数曲线，拐点对应的K就是合理的稀疏度。注意拐点不一定是全局最低点，取曲率变化最大处更稳妥。

还有一个容易踩的细节：正频率和负频率是成对出现的，设置K时不能只数正弦个数。比如三个实数余弦，频谱里实际有六条非零谱线，K最小也要给到6，否则正负频率只能恢复一半，重构出来波形完全不对。

4.2 采样点数M与恢复成功率的经验关系：M约等于2K·log(N/M)背后的直觉

理论界常写成M ≥ C·K·log(N/M)，C是某个常数。这个公式对初学者的实用价值不在精确计算，而在给出两个直觉：一是M随K线性增长，信号越稀疏越省采样；二是M随N对数增长，网格点数翻倍，采样点数不需要跟着翻倍，这也是CS最吸引人的地方。

实际调试时，我不会只算理论值，而是直接在固定N和K下扫M，观察恢复误差曲线。以N=512、K=3（三个实正弦）为例，典型行为是：M=32时偶尔能恢复，稍加噪声就翻车；M=64时稳定恢复，频率和幅度都基本正确；M=128时余量很足，换随机种子也不容易失败。

4.3 频率网格分辨率：字典不匹配时怎么降误差

DFT网格的间距是fs/N。信号真实频率不在网格上时，任何基于网格的字典算法都只能选出最近的网格点，这就是字典失配，也叫离网格问题。OMP恢复出的频率永远是网格点，真实频率和最近的网格点差多少，频率误差就有多大。

三个常用对策。一是加长数据窗，把N做大，频率网格变细，这是最朴素的做法，代价是字典矩阵变大，OMP计算量上升；二是局部细化，用粗网格先找频率的大致位置，再在这个位置附近加密网格，跑第二次OMP或直接做精细搜索；三是干脆抛弃网格，把频率当成连续变量做非线性最小二乘拟合，这就是第六章的做法。

注意一个反直觉的地方：加密频率网格不总是好事。网格越密，相邻原子的相关性越高，OMP在原子选择时越容易犹豫。当两个真实频率靠得很近时，加密网格反而加剧原子间的相干性，让恢复从“能分出两个峰”退化成“只能选一个峰”。碰到这种情况，优先考虑增加M或者改用基追踪，而不是无脑加密网格。

5. 压缩感知恢复的五个避坑记录：字典失配、直流分量与随机种子翻车

这一章专门收录我用OMP恢复正弦信号时实际碰过的典型失败案例。每条按“现象→原因→解决”展开，调试时可以直接对照。

5.1 现象：恢复出的频率整体偏移零点几个Hz——字典频率网格失配

用50.3Hz的正弦做实验，OMP恢复出的频率是50.0Hz，幅度看起来还行，但频率精确值差了0.3Hz。把信号频率改成50.0Hz，恢复又变成完美。这种现象在真实数据上非常普遍，因为真实信号的频率极少恰好落在DFT网格上。

原因是字典的频率网格只有512个原子，间隔fs/N=2Hz，真实频率落在两个网格点之间，OMP只能选内积更大的那一个，频率精度被网格限制在正负1Hz以内。这属于模型误差，不是算法bug。

解决思路分两步：先用粗网格OMP确定频率所在的大致bin区间，再这个区间内加密网格，比如把网格细化到0.01Hz，重新跑一次OMP。如果对精度要求更高，直接做第六章的连续频率精化，把频率当作连续变量求解，彻底摆脱网格限制。

5.2 现象：明明只有正弦，恢复结果却出现直流伪峰——忘了去均值

信号里有一个0.2V的直流偏置，其他部分是两个正弦。OMP恢复出的频谱在0Hz处出现一个幅度很大的峰，而且有时候这个直流峰会压制附近的低频分量，导致低频正弦恢复不出来。

原因在于直流分量在频域同样是一个稀疏原子，而且它的原子向量是常数向量，和所有观测数据的均值直接相关。只要信号均值不为零，0Hz原子就能解释掉一大块能量，OMP第一步往往就选中它。这本身没有错，但如果你以为信号是纯正弦，看到直流峰就会误判。

解决方法是先在预处理阶段去除均值，即x = x - mean(x)，然后恢复完频谱后再把均值信息单独保留。如果直流分量本身就是你要分析的目标，那就把它当作一个合法的字典原子，注意K要相应加1。观察欠采样数据时，也可以用y - mean(y)做一次粗略去偏，但严格做法是在字典里显式包含直流原子。

5.3 现象：换一次随机种子结果差很多——观测矩阵的相干性不稳定

同一个信号，rng(1)跑出来三个频率都恢复了，改成rng(2)就漏掉120Hz那一路。有时采样率M都到理论值的两倍了，换种子后依然失败。

原因是randperm每次抽取的样本集不同，观测矩阵也不同。某些采样集合会让两个频率的原子的投影过度接近，相当于观测序列里恰好缺了能区分这两个频率的关键时刻。稀疏恢复理论给的是“高概率”保证，不是“必然”保证，单次随机实例完全可能落在坏集合里。

解决思路是固定随机种子做验证；确定最终方案时，可以多试几组种子，看看哪一组恢复结果最稳。更正规的做法是设计确定性采样模式，比如用Sobol低差异序列代替纯随机抽取，让采样点在时间轴上分布得更均匀，既保留随机性又降低坏集合出现的概率。

5.4 现象：OMP在低信噪比下把噪声当信号——残差阈值设得不合理

信号加了10%幅度的高斯白噪声之后，OMP恢复出的频谱出现一堆小尖峰，分布在各个频率位置，真实谱线混在其中。这个现象的典型特征是伪峰数量随K增大一起增多，且伪峰位置每次运行都不一样。

原因是噪声在字典上的投影不可能为零，OMP每轮都在残差里找到相关性最大的原子。如果残差阈值tol设得过小，比如1e-6，算法会把噪声的投影也当作有效信号，一直迭代到K上限。

解决方法是把tol设成有意义的物理量。比如先估计噪声标准差sigma，再取tol = 3 * sigma * sqrt(M)，让触发器在残差落入噪声水平时停止。另一个思路是不设tol，改用第四章的残差拐点法，观察残差范数曲线的下降趋势，手动确定该停在哪一步。对信噪比低于10dB的场景，建议直接考虑基追踪，OMP的贪心策略在强噪声下确实吃亏。

5.5 现象：两个频率靠得近时只能恢复出一个——原子库相关性爆炸

两个正弦频率分别是100Hz和103Hz，N=512，fs=1024，理论上频率分辨率为2Hz，两个峰应该分得开。但OMP只恢复了100Hz那一个，103Hz完全消失。

原因是100Hz和103Hz对应网格原子之间的内积接近于1，两个字典原子几乎线性相关。OMP在第一轮选了100Hz原子后，新残差里103Hz的信息已经被大幅削弱，第二轮即使选中103Hz，幅度也会被严重低估，算法误以为它没有那么重要。

解决方法是先加大M，观测信息越充分，两个相似原子被区分开的可能性越大；然后考虑用基追踪，l1最优化对高度相关的字典比贪心算法更耐受；如果频率间隔小于瑞利限，那么问题已经不是压缩感知能单独解决的，需要结合子空间方法，比如用ESPRIT做局部精化。加密频率网格在这里是反效果，越密相关越强，谨记这一点。

6. 用交叉验证定稀疏度，再用最小二乘精化频率网格

最后一章给两个让恢复结果从“能用”变“可信”的手段：交叉验证定K，连续频率精化去网格偏差。这两步做完，恢复出的频率和幅度才有底气作为最终结果输出。

6.1 把采样点拆成训练集和验证集，扫K找误差拐点

不知道真实稀疏度时，最可靠的办法是把观测点拆成两份，一份参与OMP恢复，另一份单独用来验证。OMP用训练点做恢复，验证点的观测值则用来检验恢复质量，K给大了或给小了，验证误差都会上升。

idx = sort(randperm(N, M).'); % 全部采样位置 % 按 6:4 拆分训练集和验证集 Mtrain = round(M * 0.6); idx_tr = idx(1:Mtrain); idx_val = idx(Mtrain+1:end); % 分别构造观测矩阵 A_tr = exp(1i * 2 * pi * t(idx_tr) * f_grid) / sqrt(N); A_val = exp(1i * 2 * pi * t(idx_val) * f_grid) / sqrt(N); val_err = zeros(15, 1); for Ktest = 1:15 s_t = omp_recover(A_tr, y(idx_tr), Ktest, 1e-10); res = y(idx_val) - A_val * s_t; val_err(Ktest) = norm(res) / norm(y(idx_val)); end [~, bestK] = min(val_err); fprintf('交叉验证选出的稀疏度: %d\n', bestK); figure; plot(1:15, val_err, '-o'); xlabel('K 取值'); ylabel('验证集相对误差');

验证误差曲线在K过小时较高，因为模型缺少原子解释不了全部信号；K超过真实稀疏度后，误差会先下降到一个平台，再因为伪峰拟合噪声而缓慢上升。取曲线的最低点对应的K作为稀疏度，比拍脑袋给K靠谱得多。不过要注意，采样点数M本身不足时，这个曲线会变得很平坦，最低点位置也不稳定，说明当前数据支持不了更多稀疏度的猜测。

6.2 用OMP恢复出的频率做连续频率最小二乘精化

网格误差的彻底解法是把频率当作连续变量，在OMP给出的初值附近做非线性最小二乘。

function r = sin_fit_residual(fvec, t_obs, y, N) % 连续频率模型残差：频率为未知量，幅相由线性最小二乘消去 Afit = exp(1i * 2 * pi * t_obs * fvec(:).') / sqrt(N); cfit = Afit \ y; % 线性幅相系数 r = y - Afit * cfit; % 模型残差 end % 主脚本调用 f_init = f_grid(pos(1:3)); % OMP 选出的三个频率作为初值 opts = optimoptions('lsqnonlin', 'Display', 'off'); f_fit = lsqnonlin(@(f) sin_fit_residual(f, t_obs, y, N), ... f_init, f_init-0.5, f_init+0.5, opts); fprintf('精化后的频率:\n'); fprintf(' %.4f Hz\n', f_fit);

这个做法的关键在于可分离性：频率是非线性参数，但幅相是线性参数，所以每一步迭代里先用最小二乘解出幅相，残差里只保留频率相关的非线性部分。初值来自OMP，离真实解已经不远，边界只设正负0.5Hz也足够收敛，不容易陷入远处局部极值。

我现在的习惯是：任何CS恢复实验都必须跑通这两步才下结论。第一步用交叉验证确认稀疏度，第二步用连续精化把频率修准。这样做可以让输出结果有明确的物理解释——三个正弦就是三个频率、三个幅度、三个相位，而不是字典里的一堆离散原子值。这个习惯帮我避开了不少“恢复出来挺漂亮但一拿去做定量分析就出错”的尴尬场合。希望帮到你。

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