简介压缩感知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而言的。N512时稀疏度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。这正是“随机欠采样”这个名字的来源。很多初学者会问均匀间隔抽行不也是欠采样吗行不行答案是几乎肯定不行。均匀抽取会让奈奎斯特采样定理意义上的频谱混叠集中爆发——原始谱线镜像折叠到低频处和真实谱线叠在一起稀疏性直接被破坏。随机抽取则不同混叠能量像噪声一样平摊到整个频谱真实谱线仍然作为大峰值冒出来贪心算法就有机会逐步找到它们。这是随机欠采样相比均匀欠采样的本质优势。采样方式观测矩阵混叠特征对稀疏恢复的作用均匀下采样固定间隔抽行混叠集中在镜像频率上与真实谱线重叠稀疏性被破坏基本无法恢复随机欠采样randperm抽行混叠像噪声一样摊平真实谱线保持大峰值可恢复2.3 恢复算法的选型OMP、基追踪与阈值迭代的取舍稀疏恢复算法大体可分为贪心类、凸优化类和迭代阈值类。OMP是贪心类代表每次选一个与残差相关性最强的原子然后用最小二乘更新系数迭代K步。它的实现成本低跑得快在稀疏度小、信噪比高的场景下表现可靠这也是本文选择OMP作为主算法的原因。基追踪Basis Pursuit把l0问题松弛成l1凸优化对字典原子高度相关的情况更耐受但需要CVX或l1-magic这类外部优化工具部署重一点。迭代硬阈值IHT的迭代格式简单适合超大规模问题但收敛慢参数敏感性较高。三者的取舍可以简化为先跑OMP如果发现两个频率靠太近恢复不稳定再上基追踪信号维度极大时考虑IHT。算法实现难度对字典相干度敏感典型适用场景OMP低十几行Matlab较敏感低稀疏度、高信噪比、快速验证基追踪中依赖CVX耐受性好频率接近、字典网格不理想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(end1) 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未知时扫描K1到K20画残差范数曲线拐点对应的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观察恢复误差曲线。以N512、K3三个实正弦为例典型行为是M32时偶尔能恢复稍加噪声就翻车M64时稳定恢复频率和幅度都基本正确M128时余量很足换随机种子也不容易失败。采样数M采样比例实际效果N512K3326.2%临界状态成功率不可控6412.5%能恢复建议作为最低配置9618.8%稳定区域抗噪能力较好12825%余量大适合验证新算法或加噪声实验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/N2Hz真实频率落在两个网格点之间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和103HzN512fs1024理论上频率分辨率为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(Mtrain1: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_init0.5, opts); fprintf(精化后的频率:\n); fprintf( %.4f Hz\n, f_fit);这个做法的关键在于可分离性频率是非线性参数但幅相是线性参数所以每一步迭代里先用最小二乘解出幅相残差里只保留频率相关的非线性部分。初值来自OMP离真实解已经不远边界只设正负0.5Hz也足够收敛不容易陷入远处局部极值。我现在的习惯是任何CS恢复实验都必须跑通这两步才下结论。第一步用交叉验证确认稀疏度第二步用连续精化把频率修准。这样做可以让输出结果有明确的物理解释——三个正弦就是三个频率、三个幅度、三个相位而不是字典里的一堆离散原子值。这个习惯帮我避开了不少“恢复出来挺漂亮但一拿去做定量分析就出错”的尴尬场合。希望帮到你。本文还有配套的精品资源点击获取