简介这一MATLAB代码资源专注于压缩感知框架下的稀疏信号重建主要面向通信、雷达、图像处理等高年级本科生、研究生以及刚接触该方向的研发人员。资源将完整流程拆分为两部分OMP算法核心函数实现了正交匹配追踪中的原子选择、最小二乘更新和残差计算演示脚本则完成从原始稀疏信号生成、观测矩阵采样到OMP恢复再对比误差的全过程便于读者结合运行结果理解压缩感知“低维采样、高维重建”的基本思想。包内共2个文件全部为.m格式的MATLAB源码整体压缩包仅912B体量极小、无额外依赖适合在MATLAB中直接运行与逐行调试也可作为课程作业或论文复现的基础模板。目前该资源已有2247人浏览学习对于希望以低成本快速上手OMP算法并进行压缩感知基础实验的读者来说这是一份轻量且实用的参考实现。 这几年只要做信号处理方向几乎都会碰到同一个问题导师丢过来一句话“去把压缩感知做个仿真出来”然后就忙自己的事去了。我第一次拿到这个题目半天时间都耗在搞清OMP和压缩感知到底是什么关系上剩下的时间全花在让MATLAB跑出正确结果。这篇东西就是把我当时的理解、代码和后来踩过的坑一起整理出来给正打算用OMP做压缩感知信号重建的人一个能直接跑的起点。OMP全称Orthogonal Matching Pursuit正交匹配追踪是压缩感知里最容易上手的重建算法之一。配合MATLAB做一维稀疏信号重建整套流程大概半小时就能跑通。读这篇文章之前你只需要知道线性代数的基本概念我会把模型、算法步骤、完整代码和常见问题全部放在一起最后再看几个仿真结果可以对照整个过程基本没有黑盒。1. 压缩感知与OMP先弄懂为什么选这条路1.1 压缩感知到底在干什么传统信号采集遵循奈奎斯特采样定理采样率至少要达到信号最高频率的两倍否则频谱混叠信号就废了。压缩感知的思路完全不同它建立在“大部分自然信号在某个变换域下是稀疏的”这个事实上。一个时域上看似复杂的信号在傅里叶变换域或者小波变换域里可能只有少数几个非零系数其他都是零或接近零。既然信号本身携带的有效信息远少于它的长度那能不能在采集阶段就把冗余压缩掉压缩感知给出的答案是可以。它通过一个M×N的测量矩阵Φ把N维稀疏信号x直接投影到M维观测向量y上数学表达式是y Φx其中M远小于N。比如一个长度N256的信号我可能只用M64个观测值就能重建它。这在传统采样框架下是不敢想象的因为方程yΦx是一个欠定方程未知量256个方程只有64个正常来说解有无数个。但压缩感知理论告诉我们只要x足够稀疏且Φ满足约束等距性RIP就可以从这64个观测里把256个值准确重建出来。1.2 为什么入门首选OMP重建算法是压缩感知的核心。主流方法分两类一类是基追踪Basis Pursuit把问题转化为l1范数最小化凸优化典型工具有CVX、SPGL1另一类是贪婪类算法OMP就是其中最经典的一个。我自己的体会是OMP作为入门首选有三个原因。第一思路极其直观每一步都在和“相关性”打交道不像凸优化那样要理解对偶、内点法一堆东西第二MATLAB实现只要几十行不需要额外装任何工具箱对只想快速跑通验证的新手非常友好第三在小规模、严格稀疏信号场景下OMP的重建精度和速度都很能打。当然它不是万能的。基追踪在噪声较强、信号不是严格稀疏而只是可压缩的时候往往更稳代价是求解速度慢、代码复杂。OMP则在强稀疏场景下速度快、内存占用低尤其适合做实时或嵌入式方向的原型验证。对比项OMP贪婪类基追踪凸优化类实现难度低几十行代码高需要优化工具包运行速度快慢噪声鲁棒性一般较好理论保证较弱强l1最小化有理论极限适用场景严格稀疏、规模小可压缩信号、大规模、强噪声如果你只是做课程作业或者验证压缩感知的基本概念直接用OMP就够了。等后面碰到真实信号重建效果不理想再考虑换基追踪或更进阶的算法不迟。2. OMP算法核心原理与数学推导2.1 从“找原子”开始理解测量矩阵Φ可以看作一沓“原子”每一列都是一个原子一共N个。观测向量y则是这些原子在未知系数x下的线性组合。现在观测值y已知但不知道是由哪些原子组合出来的。OMP的思路就是逐一挑选“嫌疑犯”每次从这一堆原子中找出与当前残差最相关的那一列把它加入支撑集然后利用最小二乘重新计算支撑集上的系数再更新残差循环往复。“相关”在这里用内积衡量。第j个原子与残差r的内积是corr(j) φ_j, r内积绝对值越大说明这个原子越像残差的一部分。第一次迭代时残差就是y本身所以要找的是和观测向量最像的那个原子。这一点用生活化的方式理解就是你丢了一个图案手上有一堆拼图碎片第一块肯定挑轮廓最匹配的拼上去之后看剩下缺口的样子再挑下一块。OMP做的就是这个事只不过“缺口”变成了残差。2.2 OMP四步迭代流程完整算法流程可以拆成四步初始化残差r_0 y支撑集S 空集迭代次数iter 1。原子选择计算所有原子与当前残差的内积找到绝对值最大的索引idx argmax |φ_j, r_{iter-1}|把idx加入支撑集S。系数更新在支撑集张成的子空间上求解最小二乘x_s pinv(Φ(:, S)) * y注意这里用的是伪逆pinv因为Φ(:,S)是欠定矩阵普通求逆不可行。残差更新r y - Φ(:, S) * x_s然后回到第2步直到迭代次数达到稀疏度K或残差范数小于设定阈值。这里面最容易忽略的是第3步的最小二乘很多人初学时直接取x_s中对应位置等于第2步的内积值不重新计算。结果就是残差更新得不准同一个原子可能被反复选中支撑集越选越偏最后重建出来完全是乱的。2.3 为什么最小二乘这步不能省从几何上看残差更新其实是在做一次正交投影把y投影到支撑集列向量张成的子空间上残差是这个投影的补量。只有用了最小二乘求出的系数投影才是正交投影残差才会和支撑集里的每一个原子都正交这样才能保证已经选过的原子不会再被选中。我用一个比喻解释第一次选了和y最像的原子φ1如果没有最小二乘校准残差里可能还残留φ1方向的分量导致下一次又选到φ1。做了最小二乘之后φ1方向的分量被完全扣掉残差只保留其他方向的信息选择过程才有“新信息”可言。这也是OMP这个名字里“正交”两个字的由来。提示代码里如果出现支撑集索引重复大概率是残差更新时没有用最小二乘结果或者伪逆用错了对象。优先检查这一步。理解了这一点OMP的整个逻辑就通了每一步都在当前残差里找最匹配的原子找到后把已经解释过的部分完全剔除再用剩下的残差继续找。K次迭代后支撑集里就有K个位置取x在支撑集上的估计值即可。3. MATLAB完整实现与参数选择3.1 可直接运行的Demo代码有了理论基础代码就很自然了。我这里给一个完整的demo信号为长度N256的随机稀疏信号稀疏度K10观测数M64。%% 压缩感知OMP重建Demo clear; close all; clc; rng(2025); % 固定随机种子保证结果可复现 %% 1. 生成K稀疏信号 N 256; % 信号长度 K 10; % 稀疏度 x zeros(N, 1); pos randperm(N, K); % 随机选K个非零位置 x(pos) randn(K, 1); % 非零位置赋高斯随机值 %% 2. 构造测量矩阵并观测 M 64; % 观测数 Phi randn(M, N) / sqrt(M); % 高斯随机测量矩阵 y Phi * x; % 观测向量 %% 3. OMP重建 x_hat OMP(Phi, y, K); %% 4. 评估重建质量 recon_err norm(x - x_hat) / norm(x); fprintf(相对重建误差: %.3e\n, recon_err); %% 绘图 figure; subplot(2, 1, 1); stem(1:N, x, filled); title(原始稀疏信号); subplot(2, 1, 2); stem(1:N, x_hat, filled); title(OMP重建信号);OMP函数单独放在一个脚本文件里建议命名为OMP.mfunction x_hat OMP(Phi, y, K) % 正交匹配追踪重建 [~, N] size(Phi); r y; % 残差 support []; % 支撑集 x_hat zeros(N, 1); x_s []; for iter 1:K % 向量化计算所有原子与残差的内积 corr Phi * r; [~, idx] max(abs(corr)); % 防止数值误差导致重复选择 if ismember(idx, support) break; end support [support, idx]; % 最小二乘更新支撑集上的系数 x_s pinv(Phi(:, support)) * y; % 更新残差 r y - Phi(:, support) * x_s; % 残差足够小时提前停止 if norm(r) 1e-10 break; end end if ~isempty(x_s) x_hat(support) x_s; end end这里有个细节值得说计算原子与残差的内积时很多人会写for循环逐个原子求我建议直接用Phi * rMATLAB矩阵乘法一次就算出所有N个内积。N很大时这个差别非常明显N256时还感觉不到N10000时for循环会卡到你怀疑人生。3.2 参数N、M、K到底怎么定这三个参数是压缩感知仿真的地基建议一开始就把它们的关系想清楚。K是稀疏度也就是信号非零系数的个数必须提前知道或者估计。这是OMP的固有前提也是它的软肋真实场景里K往往是未知的这时要么凭经验给一个偏大值要么用残差阈值提前停止要么换稀疏度自适应的SAMP算法。仿真阶段我们直接给定K让算法在理想条件下工作把原理跑通再说。M是观测数直接决定压缩比。理论上M的下界大约是2K ln(N/K)用自然对数算一下N256、K10时2×10×ln(25.6)≈64.5所以M64刚好在这个临界值附近重建效果已经很好。如果你想更保险取M80或M100效果会进一步提升但压缩比就下来了。实际项目中M要根据噪声水平和稀疏度折中不是一个拍脑袋的数。N是信号长度代表了问题规模。N越大K不变时重建难度会略微上升因为OMP要搜索的空间变大了。同时满足理论下界的M也会变大计算量跟着涨。3.3 如果信号在时域不稀疏怎么办很多真实信号在时域上一点也不稀疏比如几个正弦叠加的波形。这时候直接把x当作稀疏信号会失败因为x本身非零元素极多。正确的做法是引入稀疏基矩阵Ψ让信号在某个变换域里稀疏x Ψθ其中θ是稀疏系数。观测模型变成y ΦΨθOMP求解的是θ得到θ_hat后再用x_hat Ψθ_hat恢复时域信号。以频域稀疏为例构造DFT稀疏基N 256; f [10, 25, 60]; % 三个频率分量 theta zeros(N, 1); theta(f 1) [1, 0.6, 0.3]; % 频域稀疏系数 Psi exp(-2i * pi * (0:N-1) * (0:N-1) / N) / sqrt(N); x Psi * theta; % 时域信号 % 观测 M 64; Phi randn(M, N) / sqrt(M); y Phi * x; % OMP重建稀疏系数 theta_hat OMP(Phi * Psi, y, 3); x_hat Psi * theta_hat;频域只有3个非零系数所以K取3。这个思路可以推广到小波基、DCT基等关键是要选一个能让信号“变稀疏”的变换域。仿真时最常踩的坑就是忘了把稀疏基乘进去直接对非稀疏信号跑OMP结果自然惨不忍睹。4. 仿真结果分析与性能验证4.1 理想无噪声情况下的重建效果直接运行上面的demo终端输出的相对重建误差通常在1e-14量级基本就是机器精度说明理想条件下OMP能实现完全重建。图上原始信号和重建信号逐点吻合稀疏位置和幅度都对得上。这个结果验证了OMP的核心逻辑K10每个非零位置都准确找到支撑集对了系数又通过最小二乘精确算出误差自然趋近于零。很多人在这一步尝到甜头之后就直接上真实数据结果一塌糊涂因为真实数据几乎不会这么听话。我建议在demo跑通之后做一个蒙特卡洛实验验证稳定性固定N256、K10在不同M取值下分别跑100次随机信号统计成功重建率。所谓成功就是相对误差小于1e-3。我自己跑过的结果大概是M30时成功率很低M50时约七成M64时接近百分百M80后稳定在百分百。这个实验能帮你直观感受观测数和重建成功率之间的关系比单次实验有用得多。4.2 添加噪声后的行为变化真实系统里观测值必然带噪声模型变成y Φx n。这时OMP的行为会发生变化如果继续固执地迭代K次它会把噪声也当作信号分量去拟合重建结果出现一堆小毛刺稀疏散点图变得乱糟糟。解决办法是把停止条件从“固定迭代K次”改成“残差范数降到阈值以下”。比如设置阈值tol 1e-3 * norm(y)当残差范数小于该阈值就直接退出。这样做的好处是噪声存在时残差下降到一定水平后进入平台期继续迭代收益很小强行迭代反而过度拟合。我实际测试过信噪比20dB左右时阈值法重建出的信号比固定迭代法平滑很多误差能小一个数量级。注意阈值也不能设得太小。设成1e-12的话带噪声时残差永远降不到这个数循环就退化成固定迭代K次设得太大会提前停止丢失幅度较小的真实分量。一般从1e-2 * norm(y)到1e-4 * norm(y)之间调试比较合适。4.3 稀疏度K对结果的影响K是OMP里最敏感的参数。K给小了比如真实稀疏度是10你只给5那只重建出5个非零位置其余信号分量全部丢失误差自然大。K给大了比如给了15算法会硬凑15个支撑位置把本地的噪声残差也解释掉产生虚假分量。所以在实际调试中如果完全不知道K我的做法是先跑一个较大的K观察残差范数随着迭代次数的下降曲线。曲线通常有一个明显的拐点拐点之前每步残差下降很快拐点之后下降速度骤减这个拐点大概率就是真实稀疏度附近。这个技巧在真实数据上非常实用算是OMP调参的一个独家经验。5. 常见问题与排查技巧实录5.1 高频问题速查表我自己和周围同学实际踩过的坑整理成了一张表现象原因排查与解决重建误差在0.5以上M过小或K估计错误增大M至2K ln(N/K)以上核对稀疏信号非零个数支撑集出现重复索引残差更新时没用最小二乘结果检查pinv求解和残差更新顺序残差下降很慢信号在时域不稀疏引入稀疏基矩阵重构优化目标为theta有噪声时重建结果毛刺多固定迭代K次把噪声一起拟合改用残差阈值停止norm(r) 1e-3 * norm(y)每次运行结果不一样测量矩阵随机生成未固定种子脚本头部加rng(2025)固定随机数生成器N较大时运行极慢用for循环计算原子内积换成Phi * r一次性计算全部相关值pinv结果异常支撑集原子近似线性相关重新生成测量矩阵或对测量矩阵做列归一化5.2 几个容易被忽略的操作细节第一个是稀疏信号生成时的维度问题。我见过不少同学写x(pos) randn(1, K)然后报维度不匹配错误。原因是randn(1,K)是行向量而x是列向量MATLAB在赋值时可能自动reshape结果非零位置的数值顺序和你预想的不一样导致稀疏位置统计错误。建议始终使用randn(K, 1)保持维度清晰。第二个是测量矩阵的缩放。默认的randn(M, N)每列范数是一个均值为sqrt(M)的随机值列与列之间范数略有差异。这个差异理论上影响内积比较的公平性但因为随机矩阵列范数相对集中通常不导致选错原子。如果你做严格对比实验最好对列做归一化处理消除这个干扰。第三个是伪逆pinv和反斜杠\的区别。pinv基于SVD分解对秩亏矩阵更稳健但速度慢反斜杠对满秩方阵或最小二乘问题更快。在OMP里由于支撑集列数通常不大于M且真实信号支撑集列近似线性无关两者结果几乎一致但pinv在出现近线性相关时不至于直接崩溃所以demo里我用pinv。追求性能时可以换成反斜杠但建议先跑通再优化。5.3 从Demo走向实际应用前的准备demo跑通只是第一步。我之前接过一个振动信号重建的活直接把demo里的理想稀疏信号换成实测数据重建结果惨不忍睹。问题不在算法而在于我把“在时域稀疏”当成了默认假设。真实振动信号在频域才有稀疏性而且不可能严格稀疏都是近似稀疏小系数一堆重建时这些“漏网之鱼”会累积成不小的误差。面对这类问题我的建议是先对信号做变换域分析画出频域或小波域系数的幅度排序图确认稀疏性到底在哪个域再用降采样的方式对重建结果做交叉验证因为原始完整信号是有的可以知道重建误差的真实大小最后再考虑用更高级的算法比如SAMP自适应估计稀疏度或者CoSaMP提升抗噪能力。OMP这套东西做下来我对压缩感知的整个流程算是真正建立了手感。它代码短、逻辑直白、每一步都能在纸上画出来作为入门第一课特别合适。如果你之后要处理真实的传感器数据记得先把稀疏基和阈值这两个环节想清楚它们才是决定重建质量的关键。我自己的习惯是每个实验都保存一份随机种子和参数记录否则同一段代码两次跑出来的结果可能差得很远排查起来会非常头疼这个细节请务必保留下来。本文还有配套的精品资源点击获取 压缩感知入门：OMP算法原理与MATLAB仿真实现
简介这一MATLAB代码资源专注于压缩感知框架下的稀疏信号重建主要面向通信、雷达、图像处理等高年级本科生、研究生以及刚接触该方向的研发人员。资源将完整流程拆分为两部分OMP算法核心函数实现了正交匹配追踪中的原子选择、最小二乘更新和残差计算演示脚本则完成从原始稀疏信号生成、观测矩阵采样到OMP恢复再对比误差的全过程便于读者结合运行结果理解压缩感知“低维采样、高维重建”的基本思想。包内共2个文件全部为.m格式的MATLAB源码整体压缩包仅912B体量极小、无额外依赖适合在MATLAB中直接运行与逐行调试也可作为课程作业或论文复现的基础模板。目前该资源已有2247人浏览学习对于希望以低成本快速上手OMP算法并进行压缩感知基础实验的读者来说这是一份轻量且实用的参考实现。 这几年只要做信号处理方向几乎都会碰到同一个问题导师丢过来一句话“去把压缩感知做个仿真出来”然后就忙自己的事去了。我第一次拿到这个题目半天时间都耗在搞清OMP和压缩感知到底是什么关系上剩下的时间全花在让MATLAB跑出正确结果。这篇东西就是把我当时的理解、代码和后来踩过的坑一起整理出来给正打算用OMP做压缩感知信号重建的人一个能直接跑的起点。OMP全称Orthogonal Matching Pursuit正交匹配追踪是压缩感知里最容易上手的重建算法之一。配合MATLAB做一维稀疏信号重建整套流程大概半小时就能跑通。读这篇文章之前你只需要知道线性代数的基本概念我会把模型、算法步骤、完整代码和常见问题全部放在一起最后再看几个仿真结果可以对照整个过程基本没有黑盒。1. 压缩感知与OMP先弄懂为什么选这条路1.1 压缩感知到底在干什么传统信号采集遵循奈奎斯特采样定理采样率至少要达到信号最高频率的两倍否则频谱混叠信号就废了。压缩感知的思路完全不同它建立在“大部分自然信号在某个变换域下是稀疏的”这个事实上。一个时域上看似复杂的信号在傅里叶变换域或者小波变换域里可能只有少数几个非零系数其他都是零或接近零。既然信号本身携带的有效信息远少于它的长度那能不能在采集阶段就把冗余压缩掉压缩感知给出的答案是可以。它通过一个M×N的测量矩阵Φ把N维稀疏信号x直接投影到M维观测向量y上数学表达式是y Φx其中M远小于N。比如一个长度N256的信号我可能只用M64个观测值就能重建它。这在传统采样框架下是不敢想象的因为方程yΦx是一个欠定方程未知量256个方程只有64个正常来说解有无数个。但压缩感知理论告诉我们只要x足够稀疏且Φ满足约束等距性RIP就可以从这64个观测里把256个值准确重建出来。1.2 为什么入门首选OMP重建算法是压缩感知的核心。主流方法分两类一类是基追踪Basis Pursuit把问题转化为l1范数最小化凸优化典型工具有CVX、SPGL1另一类是贪婪类算法OMP就是其中最经典的一个。我自己的体会是OMP作为入门首选有三个原因。第一思路极其直观每一步都在和“相关性”打交道不像凸优化那样要理解对偶、内点法一堆东西第二MATLAB实现只要几十行不需要额外装任何工具箱对只想快速跑通验证的新手非常友好第三在小规模、严格稀疏信号场景下OMP的重建精度和速度都很能打。当然它不是万能的。基追踪在噪声较强、信号不是严格稀疏而只是可压缩的时候往往更稳代价是求解速度慢、代码复杂。OMP则在强稀疏场景下速度快、内存占用低尤其适合做实时或嵌入式方向的原型验证。对比项OMP贪婪类基追踪凸优化类实现难度低几十行代码高需要优化工具包运行速度快慢噪声鲁棒性一般较好理论保证较弱强l1最小化有理论极限适用场景严格稀疏、规模小可压缩信号、大规模、强噪声如果你只是做课程作业或者验证压缩感知的基本概念直接用OMP就够了。等后面碰到真实信号重建效果不理想再考虑换基追踪或更进阶的算法不迟。2. OMP算法核心原理与数学推导2.1 从“找原子”开始理解测量矩阵Φ可以看作一沓“原子”每一列都是一个原子一共N个。观测向量y则是这些原子在未知系数x下的线性组合。现在观测值y已知但不知道是由哪些原子组合出来的。OMP的思路就是逐一挑选“嫌疑犯”每次从这一堆原子中找出与当前残差最相关的那一列把它加入支撑集然后利用最小二乘重新计算支撑集上的系数再更新残差循环往复。“相关”在这里用内积衡量。第j个原子与残差r的内积是corr(j) φ_j, r内积绝对值越大说明这个原子越像残差的一部分。第一次迭代时残差就是y本身所以要找的是和观测向量最像的那个原子。这一点用生活化的方式理解就是你丢了一个图案手上有一堆拼图碎片第一块肯定挑轮廓最匹配的拼上去之后看剩下缺口的样子再挑下一块。OMP做的就是这个事只不过“缺口”变成了残差。2.2 OMP四步迭代流程完整算法流程可以拆成四步初始化残差r_0 y支撑集S 空集迭代次数iter 1。原子选择计算所有原子与当前残差的内积找到绝对值最大的索引idx argmax |φ_j, r_{iter-1}|把idx加入支撑集S。系数更新在支撑集张成的子空间上求解最小二乘x_s pinv(Φ(:, S)) * y注意这里用的是伪逆pinv因为Φ(:,S)是欠定矩阵普通求逆不可行。残差更新r y - Φ(:, S) * x_s然后回到第2步直到迭代次数达到稀疏度K或残差范数小于设定阈值。这里面最容易忽略的是第3步的最小二乘很多人初学时直接取x_s中对应位置等于第2步的内积值不重新计算。结果就是残差更新得不准同一个原子可能被反复选中支撑集越选越偏最后重建出来完全是乱的。2.3 为什么最小二乘这步不能省从几何上看残差更新其实是在做一次正交投影把y投影到支撑集列向量张成的子空间上残差是这个投影的补量。只有用了最小二乘求出的系数投影才是正交投影残差才会和支撑集里的每一个原子都正交这样才能保证已经选过的原子不会再被选中。我用一个比喻解释第一次选了和y最像的原子φ1如果没有最小二乘校准残差里可能还残留φ1方向的分量导致下一次又选到φ1。做了最小二乘之后φ1方向的分量被完全扣掉残差只保留其他方向的信息选择过程才有“新信息”可言。这也是OMP这个名字里“正交”两个字的由来。提示代码里如果出现支撑集索引重复大概率是残差更新时没有用最小二乘结果或者伪逆用错了对象。优先检查这一步。理解了这一点OMP的整个逻辑就通了每一步都在当前残差里找最匹配的原子找到后把已经解释过的部分完全剔除再用剩下的残差继续找。K次迭代后支撑集里就有K个位置取x在支撑集上的估计值即可。3. MATLAB完整实现与参数选择3.1 可直接运行的Demo代码有了理论基础代码就很自然了。我这里给一个完整的demo信号为长度N256的随机稀疏信号稀疏度K10观测数M64。%% 压缩感知OMP重建Demo clear; close all; clc; rng(2025); % 固定随机种子保证结果可复现 %% 1. 生成K稀疏信号 N 256; % 信号长度 K 10; % 稀疏度 x zeros(N, 1); pos randperm(N, K); % 随机选K个非零位置 x(pos) randn(K, 1); % 非零位置赋高斯随机值 %% 2. 构造测量矩阵并观测 M 64; % 观测数 Phi randn(M, N) / sqrt(M); % 高斯随机测量矩阵 y Phi * x; % 观测向量 %% 3. OMP重建 x_hat OMP(Phi, y, K); %% 4. 评估重建质量 recon_err norm(x - x_hat) / norm(x); fprintf(相对重建误差: %.3e\n, recon_err); %% 绘图 figure; subplot(2, 1, 1); stem(1:N, x, filled); title(原始稀疏信号); subplot(2, 1, 2); stem(1:N, x_hat, filled); title(OMP重建信号);OMP函数单独放在一个脚本文件里建议命名为OMP.mfunction x_hat OMP(Phi, y, K) % 正交匹配追踪重建 [~, N] size(Phi); r y; % 残差 support []; % 支撑集 x_hat zeros(N, 1); x_s []; for iter 1:K % 向量化计算所有原子与残差的内积 corr Phi * r; [~, idx] max(abs(corr)); % 防止数值误差导致重复选择 if ismember(idx, support) break; end support [support, idx]; % 最小二乘更新支撑集上的系数 x_s pinv(Phi(:, support)) * y; % 更新残差 r y - Phi(:, support) * x_s; % 残差足够小时提前停止 if norm(r) 1e-10 break; end end if ~isempty(x_s) x_hat(support) x_s; end end这里有个细节值得说计算原子与残差的内积时很多人会写for循环逐个原子求我建议直接用Phi * rMATLAB矩阵乘法一次就算出所有N个内积。N很大时这个差别非常明显N256时还感觉不到N10000时for循环会卡到你怀疑人生。3.2 参数N、M、K到底怎么定这三个参数是压缩感知仿真的地基建议一开始就把它们的关系想清楚。K是稀疏度也就是信号非零系数的个数必须提前知道或者估计。这是OMP的固有前提也是它的软肋真实场景里K往往是未知的这时要么凭经验给一个偏大值要么用残差阈值提前停止要么换稀疏度自适应的SAMP算法。仿真阶段我们直接给定K让算法在理想条件下工作把原理跑通再说。M是观测数直接决定压缩比。理论上M的下界大约是2K ln(N/K)用自然对数算一下N256、K10时2×10×ln(25.6)≈64.5所以M64刚好在这个临界值附近重建效果已经很好。如果你想更保险取M80或M100效果会进一步提升但压缩比就下来了。实际项目中M要根据噪声水平和稀疏度折中不是一个拍脑袋的数。N是信号长度代表了问题规模。N越大K不变时重建难度会略微上升因为OMP要搜索的空间变大了。同时满足理论下界的M也会变大计算量跟着涨。3.3 如果信号在时域不稀疏怎么办很多真实信号在时域上一点也不稀疏比如几个正弦叠加的波形。这时候直接把x当作稀疏信号会失败因为x本身非零元素极多。正确的做法是引入稀疏基矩阵Ψ让信号在某个变换域里稀疏x Ψθ其中θ是稀疏系数。观测模型变成y ΦΨθOMP求解的是θ得到θ_hat后再用x_hat Ψθ_hat恢复时域信号。以频域稀疏为例构造DFT稀疏基N 256; f [10, 25, 60]; % 三个频率分量 theta zeros(N, 1); theta(f 1) [1, 0.6, 0.3]; % 频域稀疏系数 Psi exp(-2i * pi * (0:N-1) * (0:N-1) / N) / sqrt(N); x Psi * theta; % 时域信号 % 观测 M 64; Phi randn(M, N) / sqrt(M); y Phi * x; % OMP重建稀疏系数 theta_hat OMP(Phi * Psi, y, 3); x_hat Psi * theta_hat;频域只有3个非零系数所以K取3。这个思路可以推广到小波基、DCT基等关键是要选一个能让信号“变稀疏”的变换域。仿真时最常踩的坑就是忘了把稀疏基乘进去直接对非稀疏信号跑OMP结果自然惨不忍睹。4. 仿真结果分析与性能验证4.1 理想无噪声情况下的重建效果直接运行上面的demo终端输出的相对重建误差通常在1e-14量级基本就是机器精度说明理想条件下OMP能实现完全重建。图上原始信号和重建信号逐点吻合稀疏位置和幅度都对得上。这个结果验证了OMP的核心逻辑K10每个非零位置都准确找到支撑集对了系数又通过最小二乘精确算出误差自然趋近于零。很多人在这一步尝到甜头之后就直接上真实数据结果一塌糊涂因为真实数据几乎不会这么听话。我建议在demo跑通之后做一个蒙特卡洛实验验证稳定性固定N256、K10在不同M取值下分别跑100次随机信号统计成功重建率。所谓成功就是相对误差小于1e-3。我自己跑过的结果大概是M30时成功率很低M50时约七成M64时接近百分百M80后稳定在百分百。这个实验能帮你直观感受观测数和重建成功率之间的关系比单次实验有用得多。4.2 添加噪声后的行为变化真实系统里观测值必然带噪声模型变成y Φx n。这时OMP的行为会发生变化如果继续固执地迭代K次它会把噪声也当作信号分量去拟合重建结果出现一堆小毛刺稀疏散点图变得乱糟糟。解决办法是把停止条件从“固定迭代K次”改成“残差范数降到阈值以下”。比如设置阈值tol 1e-3 * norm(y)当残差范数小于该阈值就直接退出。这样做的好处是噪声存在时残差下降到一定水平后进入平台期继续迭代收益很小强行迭代反而过度拟合。我实际测试过信噪比20dB左右时阈值法重建出的信号比固定迭代法平滑很多误差能小一个数量级。注意阈值也不能设得太小。设成1e-12的话带噪声时残差永远降不到这个数循环就退化成固定迭代K次设得太大会提前停止丢失幅度较小的真实分量。一般从1e-2 * norm(y)到1e-4 * norm(y)之间调试比较合适。4.3 稀疏度K对结果的影响K是OMP里最敏感的参数。K给小了比如真实稀疏度是10你只给5那只重建出5个非零位置其余信号分量全部丢失误差自然大。K给大了比如给了15算法会硬凑15个支撑位置把本地的噪声残差也解释掉产生虚假分量。所以在实际调试中如果完全不知道K我的做法是先跑一个较大的K观察残差范数随着迭代次数的下降曲线。曲线通常有一个明显的拐点拐点之前每步残差下降很快拐点之后下降速度骤减这个拐点大概率就是真实稀疏度附近。这个技巧在真实数据上非常实用算是OMP调参的一个独家经验。5. 常见问题与排查技巧实录5.1 高频问题速查表我自己和周围同学实际踩过的坑整理成了一张表现象原因排查与解决重建误差在0.5以上M过小或K估计错误增大M至2K ln(N/K)以上核对稀疏信号非零个数支撑集出现重复索引残差更新时没用最小二乘结果检查pinv求解和残差更新顺序残差下降很慢信号在时域不稀疏引入稀疏基矩阵重构优化目标为theta有噪声时重建结果毛刺多固定迭代K次把噪声一起拟合改用残差阈值停止norm(r) 1e-3 * norm(y)每次运行结果不一样测量矩阵随机生成未固定种子脚本头部加rng(2025)固定随机数生成器N较大时运行极慢用for循环计算原子内积换成Phi * r一次性计算全部相关值pinv结果异常支撑集原子近似线性相关重新生成测量矩阵或对测量矩阵做列归一化5.2 几个容易被忽略的操作细节第一个是稀疏信号生成时的维度问题。我见过不少同学写x(pos) randn(1, K)然后报维度不匹配错误。原因是randn(1,K)是行向量而x是列向量MATLAB在赋值时可能自动reshape结果非零位置的数值顺序和你预想的不一样导致稀疏位置统计错误。建议始终使用randn(K, 1)保持维度清晰。第二个是测量矩阵的缩放。默认的randn(M, N)每列范数是一个均值为sqrt(M)的随机值列与列之间范数略有差异。这个差异理论上影响内积比较的公平性但因为随机矩阵列范数相对集中通常不导致选错原子。如果你做严格对比实验最好对列做归一化处理消除这个干扰。第三个是伪逆pinv和反斜杠\的区别。pinv基于SVD分解对秩亏矩阵更稳健但速度慢反斜杠对满秩方阵或最小二乘问题更快。在OMP里由于支撑集列数通常不大于M且真实信号支撑集列近似线性无关两者结果几乎一致但pinv在出现近线性相关时不至于直接崩溃所以demo里我用pinv。追求性能时可以换成反斜杠但建议先跑通再优化。5.3 从Demo走向实际应用前的准备demo跑通只是第一步。我之前接过一个振动信号重建的活直接把demo里的理想稀疏信号换成实测数据重建结果惨不忍睹。问题不在算法而在于我把“在时域稀疏”当成了默认假设。真实振动信号在频域才有稀疏性而且不可能严格稀疏都是近似稀疏小系数一堆重建时这些“漏网之鱼”会累积成不小的误差。面对这类问题我的建议是先对信号做变换域分析画出频域或小波域系数的幅度排序图确认稀疏性到底在哪个域再用降采样的方式对重建结果做交叉验证因为原始完整信号是有的可以知道重建误差的真实大小最后再考虑用更高级的算法比如SAMP自适应估计稀疏度或者CoSaMP提升抗噪能力。OMP这套东西做下来我对压缩感知的整个流程算是真正建立了手感。它代码短、逻辑直白、每一步都能在纸上画出来作为入门第一课特别合适。如果你之后要处理真实的传感器数据记得先把稀疏基和阈值这两个环节想清楚它们才是决定重建质量的关键。我自己的习惯是每个实验都保存一份随机种子和参数记录否则同一段代码两次跑出来的结果可能差得很远排查起来会非常头疼这个细节请务必保留下来。本文还有配套的精品资源点击获取