简介压缩感知Compressed Sensing, CS为突破奈奎斯特采样限制提供了全新思路。面向研究随机欠采样与稀疏重构的信号处理学者和工程师这份MATLAB源码包内置两个演示程序分别调用正交匹配追踪OMP和SPGL1两种算法覆盖从随机欠采样到信号恢复的完整流程便于对比不同算法在噪声场景下的重构效果与抗噪性能。包内共36个文件以25个m脚本为主涵盖算法主程序、参数配置与示例调用辅以配套的c/h源文件、编译好的mex动态库支持Windows/Linux/macOS以及readme、changelog等说明文档整体仅67KB轻量便携。已有1717人浏览学习。通过运行示例读者能直观理解随机采样模式、稀疏基选择与恢复质量的关系掌握OMP和SPGL1的实际调用方法是压缩感知入门与算法验证的实用工具适合用于课程设计、算法对比研究或工程预研。1. 压缩感知为何随机欠采样还能恢复多个正弦信号在工程里我们习惯了奈奎斯特采样定理采样频率必须高于信号最高频率的两倍否则频谱混叠。但压缩感知给了一个反直觉的结论——只要信号在某个变换基下足够稀疏用低于奈奎斯特率的随机采样点配合 L1 范数最优化也能精确恢复原信号。这并不违背信息论因为随机采样把混叠变成了近似不相干的噪声而稀疏性让信号能从噪声里“跳”出来。这个压缩包正是把这件事跑给 MATLAB 看生成多个正弦信号随机欠采样到一个比奈奎斯特低得多的采样率再通过 OMP 与 SPGL1 两种恢复算法把原始频谱找回来。适合正在做压缩感知课程作业、算法对比或者想快速看到随机欠采样到底怎么恢复信号的人。2. 压缩感知三大支柱稀疏基、测量矩阵与恢复条件2.1 稀疏表示正弦信号在频域的稀疏性压缩感知成立的第一条是信号在某个变换域里系数个数远小于信号长度。多个正弦信号叠加在频域只有少数几根谱线其余位置都是零这类信号天然满足稀疏性要求。N 256; % 信号长度 K 5; % 正弦个数 s_true zeros(N, 1); % 频域系数向量只有K个非零 f randi(N, 1, K) - 1; % 随机频率下标范围0~255 amp randi([1, 10], 1, K); % 随机幅值 phase 2 * pi * rand(1, K); % 随机初相 for i 1:K s_true(f(i)1) amp(i) * exp(1i * phase(i)); end x ifft(s_true); % 反变换到时域得到多正弦叠加信号这里s_true是频域稀疏表示x是时域信号。若按奈奎斯特思路N 个点一个都不能少而压缩感知的目标是用远小于 N 的 M 个测量点恢复s_true再通过 IFFT 得到x。注意用ifft而不是fft这样x Psi * s_true其中Psi是 IFFT 矩阵后面构造观测矩阵时保持同一套归一化约定可以少踩一个缩放坑。2.2 随机欠采样观测矩阵的构造与约束等距性第二条是观测方式。压缩感知不要求均匀采样而是随机抽取一部分时域样本。上面代码里的x是完整信号实际采集只拿到 M 个随机位置的值。随机下标取出这些值后再结合傅里叶原子构造观测矩阵 A使得y A * s_true。M 80; % 随机采样点数采样率约0.31 M_idx randperm(N, M); % 从N个位置中随机选M个 M_idx sort(M_idx); % 按时间顺序排列方便后续画图 y x(M_idx); % 实际观测值 % 构造部分傅里叶观测矩阵A每一行对应一个采样时刻的傅里叶原子 A exp(1i * 2 * pi * (0:N-1) .* M_idx / N) / sqrt(N);参数上M 决定了采样率 M/N。M 越大恢复越容易但压缩比下降M 太小则观测方程欠定程度严重恢复成败开始变得随机。RIP 条件给出的经验边界是 M 通常取稀疏度 K 的 35 倍再乘log(N/K)这个包里的 N256、K5M80 对这个量级是安全的。随机欠采样的核心在于随机测量矩阵与固定的傅里叶稀疏基之间有很高的概率满足不相干性这正是压缩感知能成立的测量层面的原因。2.3 恢复算法选型OMP 与 SPGL1 为什么能互补得到 y 和 A 之后问题变成从欠定方程中恢复 s。OMP 是贪心算法每次迭代挑选与当前残差最相关的一列加入支撑集后做最小二乘更新迭代 K 次即可。SPGL1 则把问题写成基追踪或 LASSO 形式用谱投影梯度法解凸优化不依赖稀疏度 K 的精确值抗噪性更好。对比项OMPSPGL1原理贪心迭代谱投影梯度是否依赖稀疏度 K是否抗噪性弱强计算代价低高对测量矩阵条件数敏感度高低多正弦信号场景里如果两个频率靠得很近OMP 容易因为 A 的两列相关性高而选错支撑集SPGL1 通过全局优化规避了这个问题但代价是迭代时间更长。压缩包同时提供两种算法目的就是让你在同一个信号上直观看到贪心与凸优化两条路线的差异。3. 代码包解析CS_OMP.m 与 CS_SPGL1.m 内部逻辑3.1 先认识压缩包里的 SPGL1 工具包解压后你会看到 spgl1_1_7 文件夹这是 SPGL1 v1.7 的完整发行版。里面既有 spgl1.m 主函数也有 spg_bpdn.m、spg_lasso.m、spg_bp.m、spg_mmv.m 这些上层接口以及 NormL1_project.m、NormL12_project.m 等投影算子。各文件分工大致如下文件名作用spgl1.m主入口求解一般 L1 最小化问题spg_bpdn.m基追踪去噪BPDN封装spg_lasso.mLASSO 问题封装spg_bp.m标准基追踪BP封装spg_mmv.m多测量向量MMV问题接口在 MATLAB 里使用前先把路径加进去addpath(genpath(spgl1_1_7)); savepath;addpath(genpath(...))会递归添加子目录避免漏掉 private 文件夹里的辅助函数。savepath 是否成功取决于对 MATLAB 主路径的写权限如果提示权限问题可以改用设置路径对话框保存不保存也不会影响当前会话。这一步不做后面运行 CS_SPGL1.m 时 MATLAB 会直接报“未定义函数或变量 spgl1”。3.2 CS_OMP.m正交匹配追踪的核心迭代从文件名看CS_OMP.m 应该是打包者自己封装的独立函数不依赖额外工具箱。OMP 的标准实现核心就三步找相关列、最小二乘、更新残差。function x_hat CS_OMP(y, A, K) [M, N] size(A); x_hat zeros(N, 1); r y; % 残差初始化为观测值 T []; % 支撑集索引初始为空 x_temp zeros(N, 1); for iter 1:K proj abs(A * r); % 每一列与残差的相关性 [~, idx] max(proj); % 找相关性最大的原子索引 T union(T, idx); % 加入支撑集 x_temp(T) A(:, T) \ y; % 在支撑集上做最小二乘 r y - A * x_temp; % 更新残差 end x_hat x_temp; end这个实现里最需要关注的参数是 K即稀疏度。如果 K 比真实稀疏度大OMP 会把本应是零的系数强行拟合出来引入伪峰如果 K 比真实值小残差降不下去主频率恢复不全。如果你不确定 K常见做法是先跑一次 SPGL1 看输出支撑集大小再用这个值回传给 OMP这样两种算法能互相校核。3.3 CS_SPGL1.mBPDN/LASSO 的封装逻辑SPGL1 官方主函数签名是[x, r, g, info] spgl1(A, b, tau, sigma, options);其中 tau 是 L1 范数上界sigma 是噪声能量上界。只给 tau 就是基追踪去噪BPDN只给 sigma 就是 LASSO。CS_SPGL1.m 大概率是把 y、A 和一个误差容限参数映射到 spgl1 入口上一个典型的封装思路如下function x_hat CS_SPGL1(y, A, tol) options spgSetParms(verbosity, 0); % 关闭迭代打印 [x_hat, ~, ~, info] spgl1(A, y, 0, tol, options); end这里的 tol 是残差容限不是迭代步长。如果把 tol 设为 0会变成标准基追踪对含噪信号会把噪声也拟合进去所以工程上通常设成略大于噪声能量的小数比如 1e-3。SPGL1 对 A 的列缩放很敏感A 每一列范数应当接近否则谱投影梯度的收敛会非常慢这是使用各种 L1 求解器时最常见的隐形坑。4. 完整复现生成多正弦信号、随机欠采样与恢复对比4.1 生成信号、随机欠采样、构造观测矩阵把前面的理论连成一条可跑通的脚本。先用固定随机种子保证结果可复现rng(42); N 256; % 信号长度 K 5; % 稀疏度正弦个数 M 80; % 随机采样点数 f randi(N, 1, K) - 1; % 频率下标范围0~255 amp randi([1, 10], 1, K); % 幅值1~10 phase 2*pi*rand(1, K); % 随机初相 s_true zeros(N, 1); for i 1:K s_true(f(i)1) amp(i) * exp(1i*phase(i)); end x_true ifft(s_true); % 随机欠采样从N个时域位置里选M个 M_idx randperm(N, M); M_idx sort(M_idx); y x_true(M_idx); % 构造观测矩阵A满足 y A * s_true A exp(1i * 2 * pi * (0:N-1) .* M_idx / N) / sqrt(N); % 验证模型一致性理论值应接近机器精度 disp(norm(A * s_true - y));最后一行验证非常关键。如果这个差值不是 1e-12 量级说明 A 的构造有问题后面所有恢复结果都不可信。我一般会在每次修改采样率后都跑一遍这行验证避免把时间浪费在恢复算法的错误排查上。4.2 调用 OMP 与 SPGL1 恢复数据接下来直接调用压缩包里的两个入口函数。如果你的包内函数签名与本例不同按实际接口调整参数即可s_omp CS_OMP(y, A, K); tol 1e-3; s_spgl CS_SPGL1(y, A, tol); % 反变换回时域 x_omp ifft(s_omp); x_spgl ifft(s_spgl); % 相对误差评估 err_omp norm(x_true - x_omp) / norm(x_true); err_spgl norm(x_true - x_spgl) / norm(x_true); fprintf(OMP重构误差: %.3e\n, err_omp); fprintf(SPGL1重构误差: %.3e\n, err_spgl);调用前先执行which CS_OMP确认 MATLAB 找到的是压缩包内的文件避免被工作区里其他同名旧脚本污染。相对误差能说明整体恢复精度但不能展示结构差异所以还建议画频谱对比。4.3 怎么判断恢复结果是否可信误差小不等于恢复成功还要看支撑集是否命中真实频率。画频谱叠合图stem(0:N-1, abs(s_true), b); hold on; stem(0:N-1, abs(s_omp), r); hold off; legend(真实频谱, OMP恢复);如果 OMP 的谱尖刺位置与真实谱完全重合说明稀疏位置选对了如果恢复谱里出现大量小伪峰说明 M 不够或者 A 的相干性偏高。SPGL1 恢复的谱通常更干净但 L1 软阈值特性会在非零频率附近产生很低的底噪这不影响主频率判断别因为底噪误以为恢复失败。需要特别说明的是本示例用复指数表示正弦频域只有正频率一根线逻辑最清晰。如果换成实正弦信号频域会正负对称真实稀疏度会翻倍OMP 的 K 参数也要跟着调整否则一定会漏频率。5. 压缩感知MATLAB实践避坑五个真实翻车现场5.1 现象CS_SPGL1 报“未定义函数或变量 spgl1”这是我拿到压缩包后遇到的第一件事。原因很简单spgl1_1_7 文件夹没有加入 MATLAB 搜索路径或者当前工作目录不在解压目录下。解决方法是先运行addpath(genpath(spgl1_1_7))再用which spgl1确认返回路径确实指向包内文件。不要图省事把 spgl1.m 直接复制到其他目录因为它的 private 子文件夹里还有依赖函数拆开会丢功能。5.2 现象OMP 恢复结果全是 NaN 或 Inf现象是s_omp里有 NaN画出来的频谱没法看。原因多半是 A 的条件数太大A(:, T) \ y在支撑集接近奇异时得不到稳定解尤其采样点少或两个频率靠得很近时容易发生。解决方法是把最小二乘换成 QR 分解或用伪逆pinv(A(:,T)) * y。我一般会写x_temp(T) pinv(A(:, T)) * y;pinv比\慢一些但稳定性好很多。多正弦稀疏度低支撑集尺寸通常不超过十几个这点计算代价可以接受。5.3 现象SPGL1 迭代几百步不收敛重构结果依然发散有一次我把 tol 设成 1e-8结果 SPGL1 跑了五百步还在磨蹭。原因是 tol 太小已经逼近浮点精度算法在数值噪声里原地打转。另一个常见原因是 A 的列没有归一化。解决方法是先做列归一化A A ./ vecnorm(A);再把 tol 设到噪声能量量级比如 1e-2 到 1e-4。用options spgSetParms(verbosity, 1)打开迭代打印观察前后两次目标函数值变化如果连续几十步变化很小就该调 tol 而不是继续加大迭代次数。5.4 现象换一个随机种子恢复结果天差地别随机欠采样的不同采样位置会带来不同的恢复质量这是压缩感知理论的固有概率属性不是代码 bug。一次跑通不能说明算法可靠一次失败也不能说明算法不行。正确做法是固定rng做单次调试然后在评估阶段跑 50 次蒙特卡洛统计恢复误差的中位数和失败率。我在对比 OMP 和 SPGL1 时都是分别跑 50 个种子然后画误差箱线图这样才能看出两个算法的真实差距。5.5 现象恢复出的时域波形形状对但幅值整体偏小这是缩放关系不一致的典型症状。A 的构造漏了1/sqrt(N)或者 s_true 与 x_true 之间用了fft而观测矩阵按ifft约定构造都会导致模型关系差一个常数因子。解决方法是写一行验证代码assert(norm(A * s_true - y) 1e-10);这行断言应该嵌在你的主脚本里。从那以后我再也没被这种缩放问题暗算过每次构造完观测矩阵都会强制跑一遍。6. 进阶从多正弦到实信号如何验证恢复质量与调参验证恢复质量时除了相对误差我更看重支撑集匹配率——恢复出的非零位置与真实位置有多少重合。在稀疏度已知的情况下可以这样算[~, idx_true] sort(abs(s_true), descend); [~, idx_omp] sort(abs(s_omp), descend); overlap length(intersect(idx_true(1:K), idx_omp(1:K))) / K;这个指标比相对误差更贴近压缩感知的本质稀疏支撑集找对了幅值误差可以靠后续最小二乘修正支撑集找错了误差再小也是假的。如果你想把这套代码往下游延伸可以把复指数信号换成实正弦或者任意在傅里叶基下稀疏的信号核心改动只有稀疏表示和 K 的取值。调参方面我最常用的是扫参法M 从 40 扫到 120每个 M 跑 20 次随机种子记录平均相对误差和支撑集匹配率。这样你能画出一条“采样点数量-恢复成功率”曲线不同算法的最小安全采样数一眼就能看出来。对 OMP 来说K 给大一点也许还能容忍对 SPGL1 来说tol 与噪声水平的关系才是决定恢复质量的关键。我自己的习惯是拿到这类 CS 代码先不改算法直接跑通随机欠采样基线然后把 M、K、噪声三个变量逐一扫一遍。有一回在某个稀疏成像项目里我图省事没有归一化 ASPGL1 跑出来的图像细节全是错的折腾了一整天才发现是观测矩阵缩放问题。从那以后我每次建完观测矩阵都要先验证norm(A*s_true - y)在 1e-10 量级再进入恢复环节。希望帮到你。本文还有配套的精品资源点击获取