简介：面向需要掌握稀疏表示与压缩感知的MATLAB学习者与研究人员，代码包完整实现了正交匹配追踪（OMP）算法。资源共4个文件，其中3个.m脚本分别承担算法主函数、调用示例与辅助演示，另有1张流程图或效果图辅助理解迭代选原子与残差更新过程，压缩包仅5KB，轻量紧凑，方便直接阅读和调试。目前已有272人学习/下载，特别适合信号去噪、图像恢复、特征选择等场景的初学者快速上手。通过运行示例代码，读者可以清晰看到字典构造、内积匹配、正交投影和停止判断等关键步骤，并据此扩展至MRI重建、音频处理等实际问题。整体逻辑清晰，注释简明，便于学习者根据自身需求进行改造与二次开发。
1. 自己在matlab里实现omp，为什么比下载现成工具包更稳？
压缩感知里，Orthogonal Matching Pursuit（OMP）经常被当作第一个上手的稀疏恢复算法。很多人的做法是从代码分享平台下载一个 omp_matlab 函数，然后直接调用 [x,supp] = omp(A,y,K)。这个流程在无噪声、信号严格 K 稀疏时很容易跑通；可一旦把测量矩阵换成不相干性较差的过完备字典，或者观测向量带了噪声，不同版本的实现给出的结果就开始分叉。有的版本把 K 当作真实稀疏度，有的把它当作最大迭代次数；有的会对残差阈值自动处理，有的不会。自己维护一个最小 omp_matlab 实现虽然只是二十行循环，但每个参数的行为都可控，而且把算法改写成 OLS、CoSaMP，或者接入字典学习流程时也更顺手。这篇文章适合正在做压缩感知、稀疏信道估计、图像稀疏表示，或者想在 matlab 里快速验证贪婪类算法的工程人员。
2. 从y = Ax到支撑集：OMP的迭代逻辑与停止条件
2.1 稀疏恢复模型与测量矩阵构造
处理的线性观测模型是 y = A x + n，其中 A 是 m×n 测量矩阵，实际场景里 m 经常明显小于 n，x 的非零位置只有 K 个，这就是 K 稀疏的含义。问题是当 m < n 时，A x = y 存在无穷多个解，不引入稀疏约束就无法确定地求解。压缩感知理论解决的是“在什么样的 A 和多大的 K 下，这个解可以被稳定恢复”，而 OMP 是在满足类似条件时，用贪心策略去逼近这个解。
在 matlab 里构造合成数据，第一步建议把测量矩阵列归一化。原因是 OMP 每轮比较的是“哪个原子方向与当前残差更接近”，如果某列的模比其他列大，内积自然偏大，选择结果就会偏向大范数列。列归一化的代码很简洁：
m = 128; n = 256; K = 10; A = randn(m, n); A = A ./ sqrt(sum(A.^2, 1)); % 每列向量的2-范数归一为1 x_true = zeros(n, 1); x_true(randperm(n, K)) = randn(K, 1); y = A * x_true;randn 生成的是高斯随机测量矩阵，在压缩感知验证里很常用，因为它满足有限等距性质的概率高；sum(A.^2,1)对每一列求平方和，再开方得到 1×n 的列范数向量，逐列广播后就把每列归一成单位列向量。randperm(n,K)随机指定 K 个非零位置，恢复算法的任务就是把这些位置找回来。
2.2 迭代循环：原子选择、最小二乘更新与残差正交化
第 t 次迭代开始前，支撑集是 S_t-1，残差是 r_t-1。算法先在字典所有列上做一次内积，用 matlab 写就是两行：
corr = A' * r; [~, idx] = max(abs(corr));A' * r 得到 n 个内积值，绝对值大小代表每一列对残差的解释能力。选最大值对应的索引 idx，把 A 的第 idx 列加入支撑集。
加入新列之后，需要重新计算整个支撑集上的系数，而不是只更新新加的那一个系数。因为新增一列之后，先前的最小二乘解已经不是整个支撑集上的最优解。这一步在 matlab 里用反斜杠运算符实现：
x_LS = A(:, supp) \ y; r = y - A(:, supp) * x_LS;反斜杠在 m > S 且列满秩时给出的是最小二乘解，等价于对 y 向支撑集张成的子空间做正交投影。残差 r 满足 A_supp' * r = 0，也就是说已选原子与当前残差的内积严格为零。这就是“正交匹配追踪”里“正交”二字的来源：下一次迭代再做 max(abs(A'*r)) 时，支撑集里的索引不会再以明显优势被选中，算法不会反复选同一列。对比基本匹配追踪（MP），MP 只做原子方向的减法而不做正交化，同一原子可能被重复选中，收敛路径明显曲折。
2.3 停止条件在matlab代码里的落地方式
OMP 有三种常见的停止方式，实际工程里经常是组合使用：
| 停止条件 | matlab写法 | 适合场景 | 容易踩的坑 |
|---|---|---|---|
| 固定迭代 K 次 | for 循环到 K 即停 | 已知信号严格 K 稀疏 | 噪声场景下 K 给大，噪声被当信号选进来 |
| 绝对残差阈值 | norm(r) < eps | 无噪合成实验 | eps 给太小等于没设，给太大又提前停 |
| 相对残差阈值 | norm(r) <= tol*norm(y) | 绝大多数工程场景 | tol 需随信噪比调整，不可一套参数通吃 |
下文给出的最小实现采用“K 上限 + 相对残差”双条件。原因是实际工作中很难提前知道精确稀疏度，所以 K 先给一个偏大的上界，让相对阈值在支撑集进入噪声区间时提前截断。tol 不要给得太小，否则阈值条件形同虚设，循环一定会跑满 K 次。
3. 手写一个最小可运行的 omp_matlab 核心循环
3.1 完整函数代码
以下是我在项目里惯用的最小实现，不依赖任何工具箱，复制到 .m 文件即可运行：
function [x_hat, supp] = omp_matlab(A, y, K, tol) % OMP 正交匹配追踪最小实现 % 输入: % A : m*n 测量矩阵，建议先做列归一化 % y : m*1 观测向量 % K : 最大迭代次数 / 稀疏度估计上界 % tol : 相对残差停机阈值，默认 1e-6 % 输出: % x_hat : n*1 稀疏恢复信号 % supp : 被选中原子的索引列向量 if nargin < 4 || isempty(tol) tol = 1e-6; end [m, n] = size(A); r = y; % 残差初始化为观测向量 supp = []; % 支撑集索引，初始为空 x_hat = zeros(n, 1); % 预置输出向量 for iter = 1 : K corr = A' * r; % 当前残差在字典每个原子上的投影 [~, idx] = max(abs(corr)); % 取绝对值最大对应的原子索引 if any(supp == idx) % 已选列再次被选中，说明残差已正交化 break; end supp = [supp; idx]; % 扩展支撑集 x_LS = A(:, supp) \ y; % 支撑集上的最小二乘系数，反斜杠自动走QR路径 r = y - A(:, supp) * x_LS; % 更新残差，使残差与所有已选列正交 if norm(r) <= tol * norm(y) % 相对残差达到阈值则提前退出 break; end end x_hat(supp) = x_LS; % 只在支撑集位置写回系数 end这段代码的关键点有三个。第一，初始残差直接取 y，第一轮选中的必然是观测向量在字典上投影最大的原子；第二，每轮支撑集扩张后都要完整重算最小二乘，而不是只更新新增原子方向的系数；第三，用相对残差 norm(r)/norm(y) 作为停机判断，避免了 y 本身幅值变化对阈值的影响。代码里加的重复索引检查是为了应对原子高度相干时的边界情况，理论上如果残差严格正交，同一列不会被再次选中，但浮点误差可能带来极小的非零投影。
3.2 输入参数与常见取值边界
| 参数 | 维数 | 常见取值 | 说明 |
|---|---|---|---|
| A | m×n double | m ≈ 4K~8K | 列必须归一化，否则内积选择会被大范数列带偏 |
| y | m×1 double | 均值先归零 | 含直流分量时先做去均值，避免第一轮选错原子 |
| K | 正整数 | 真实稀疏度的 1.5~2 倍 | 上界给大一点，配合 tol 截断比给小更稳妥 |
| tol | double 标量 | 无噪声 1e-10，有噪声 1e-2~1e-3 | 噪声越大 tol 越大，过小会把噪声当成信号 |
实际调参时最常犯的错误是把 K 精确设成已知的稀疏度。在无噪声理想条件下这没问题，但工程数据很少严格稀疏，K 稍大一点反而给了残差阈值发挥作用的空间。与之相对，tol 也不能设成 1e-10 想当然，含噪场景下这个阈值几乎不可能被满足，循环会一直跑到 K 次，把噪声原子也收进支撑集。
3.3 用模拟数据验证恢复正确性
验证代码可以这样写：
rng(42); m = 128; n = 256; K = 10; A = randn(m, n); A = A ./ sqrt(sum(A.^2, 1)); % 列归一化 x_true = zeros(n, 1); pos = randperm(n, K); x_true(pos) = 3 * randn(K, 1); % 放大系数幅值，便于观察恢复误差 y = A * x_true; [x_hat, supp] = omp_matlab(A, y, K + 5, 1e-10); rel_err = norm(x_hat - x_true) / norm(x_true); recall = length(intersect(supp, pos)) / K; fprintf('相对误差: %e, 支撑集重合率: %.2f\n', rel_err, recall);调用时把 K 上限设置为 K+5，是为了展示“上界 + 阈值”双条件的作用：多出来的 5 次迭代允许 OMP 继续判断是否有新原子可加，但 1e-10 的相对阈值保证了残差已经足够小时不会继续无聊迭代。正常情况下输出支撑集重合率 1.00，相对误差在 1e-12 量级，差别取决于反斜杠求解的数值精度。如果重合率低于 1.00，优先检查 A 是否列归一化，以及 y 是否带了未预料的直流分量。
4. 噪声环境下调优 omp_matlab：阈值、Gram 预计算与原子相干性
4.1 噪声下停止阈值按噪声方差估计
带噪声的观测模型是 y = A x + n，其中 n 的元素服从 N(0, σ²)。当支撑集选得正确时，残差里主要剩下的就是测量噪声，因此残差范数理论值在 σ·sqrt(m) 附近。如果 tol 设得比这个水平小很多，OMP 就会把噪声方向的投影当成“值得选择”的原子，继续扩张支撑集，最终在真实支撑之外多选几个错误位置。反过来，tol 设得太大又会在真实原子还未全部进入支撑集时提前收手。
一个可靠的阈值估计是：
sigma_n = 0.05; % 观测噪声标准差，可从系统前端估计 y_noisy = y + sigma_n * randn(m, 1); tol_noise = sqrt(m) * sigma_n / norm(y_noisy); [x_hat_noisy, supp_noisy] = omp_matlab(A, y_noisy, 2*K, tol_noise);sqrt(m) * sigma_n 是残差能量的量级估计，除以 norm(y_noisy) 是为了和函数内部的相对残差判断对齐。如果 σ_n 无法直接获得，可以用一个工程上常用的稳健近似：
sigma_hat = median(abs(A' * y)) / 0.6745;0.6745 来自正态分布的 median 与标准差的比例关系，A' * y 中那些与真实支撑集无关的投影分量近似服从零均值高斯分布，所以 median 能剥离掉大系数原子的影响。这个估计比手工试 tol 稳定得多。
4.2 用预计算 Gram 矩阵缩短 matlab 循环耗时
在 m、n 都较大的仿真里，OMP 每次迭代都要算一次 A' * r，这是一次 m×n 的矩阵乘向量的运算，总体耗时大约等于 n·m·K 次浮点运算。如果字典 A 固定不变，可以预计算 Gram 矩阵 G = A' * A，并用迭代关系跳过 A' * r：
function [x_hat, supp] = omp_matlab_gram(A, y, K, tol) [m, n] = size(A); Aty = A' * y; G = A' * A; % 预计算 Gram 矩阵，n 较大时注意内存占用 r = y; supp = []; for iter = 1:K if iter == 1 corr = Aty; % 第一轮直接用 A'*y else % r = y - A(:,supp)*x_LS，所以 A'*r = Aty - G(:,supp)*x_LS corr = Aty - G(:, supp) * x_LS; end [~, idx] = max(abs(corr)); supp = [supp; idx]; x_LS = A(:, supp) \ y; r = y - A(:, supp) * x_LS; if norm(r) <= tol * norm(y) break; end end x_hat = zeros(n, 1); x_hat(supp) = x_LS; end推导逻辑很简单：残差 r 始终等于 y - A(:,supp)·x_LS，因此 A'·r = A'·y - A'·A(:,supp)·x_LS = Aty - G(:,supp)·x_LS。每一轮只需一次 n×|S| 的矩阵乘，而 |S| 远小于 m，所以在 K 远小于 m 时收益明显。G 的存储代价是 O(n²)，n=4096 时约 134MB，n=8192 时超过 500MB，内存不足时不要硬上，这也是我把这个版本独立成函数而不是直接改写核心函数的原因之一。
4.3 原子相干性与支撑集误判的排查
OMP 是贪心算法，它的失败模式非常集中：字典里两个原子高度相关时，算法选了其中一个，残差投影会让另一个也表现得很强，很容易把错误位置提前收入支撑集。字典相干系数定义为所有互异列内积绝对值中的最大值，在 matlab 里一行代码就能检查：
G0 = abs(A' * A); G0(1:n+1:end) = 0; % 对角元是原子自身的内积，置零排除 mu = max(G0(:)); fprintf('字典相干系数 mu = %.3f\n', mu);当 mu 超过 0.9 时，OMP 选错原子的概率显著上升，尤其在信号非零系数幅值有差异的场景。实际例子是图像稀疏表示里常用的过完备 DCT 字典：列数比维度大几倍，相邻频带的原子相干性很容易接近 1，OMP 会把支撑集选在相邻的冗余基上，恢复结果看着误差大但每次都稳定得诡异。matlab 图像处理的这类任务里，我一般会先用 OMP 跑一版看支撑集索引是否频繁跳动，如果跳动明显，就改走 OLS 或者换用稀疏贝叶斯类算法。matlab 优化工具箱里的 lsqlin 也能作为子问题求解器替换反斜杠，但收敛性质并不会因此变好，真正要解决的是字典设计，不是求解器。
5. 用部分傅里叶矩阵验证omp恢复：从观测到支撑集重合率
5.1 构造部分傅里叶观测矩阵
部分傅里叶矩阵是压缩感知里非常典型的一类测量矩阵，常用于 MRI 重建和雷达回波稀疏恢复。它从 n 点 DFT 矩阵中随机抽取 m 行，对应物理含义是“只观测部分频谱”。构造代码：
rng(7); n = 1024; x_true = zeros(n, 1); x_true([120 300 450]) = [1.5 -2.0 1.0]; % 三个稀疏脉冲 sel = sort(randperm(n, 256)); % 随机选 256 个频点作为观测 if exist('dftmtx', 'file') % Signal Processing Toolbox F = dftmtx(n) / sqrt(n); else F = exp(-2j*pi*(0:n-1)'*(0:n-1)/n) / sqrt(n); end A = F(sel, :); % 部分傅里叶测量矩阵 y = A * x_true;这里 A 的行数只有 256，列数是 1024，属于明确的欠定系统。DFT 矩阵除以 sqrt(n) 后每列都是单位范数，所以不需要再做列归一化。sel 排序的意义在于让观测矩阵的行顺序不影响 matlab 的稀疏运算缓存，不排序也可以。
5.2 用 omp_matlab 做恢复并计算支撑集重合率
恢复与评估代码：
K_true = length(find(abs(x_true) > 1e-6)); % 真实支撑集大小 [x_hat, supp] = omp_matlab(A, y, K_true, 1e-10); true_pos = find(abs(x_true) > 1e-6); hit = length(intersect(supp, true_pos)) / K_true; snr_rec = 20 * log10(norm(x_true) / norm(x_hat - x_true)); fprintf('支撑集重合率: %.0f%%, 恢复信噪比: %.2f dB\n', hit*100, snr_rec);在无噪声条件下，OMP 能在这类部分傅里叶矩阵上精确恢复三个脉冲，支撑集重合率应为 100%。如果出现漏选，先确认观测频点数量是否过低，m 明显小于 4K·log(n/K) 时恢复失败属于正常情况，不是代码错误。
5.3 含噪场景下用残差能量检验支撑集可信度
最后一个实用技巧：在含噪环境里拿到支撑集之后，不要只看恢复信噪比，要单独检验支撑集是否可信。做法是保留前文估计出的 σ_n，然后比较最终残差能量与噪声理论能量：
sigma_n = 0.02; y_noisy = y + sigma_n * randn(size(y)); [x_hat, supp] = omp_matlab(A, y_noisy, K_true + 5, 1e-8); x_check = zeros(n, 1); x_check(supp) = A(:, supp) \ y_noisy; r_final = y_noisy - A(:, supp) * x_check(supp); res_energy = norm(r_final)^2; noise_energy = sigma_n^2 * length(y_noisy); if res_energy > 2 * noise_energy disp('残差偏大，支撑集可能有漏选，需要放开K或放宽tol'); elseif res_energy < 0.1 * noise_energy disp('残差异常小，存在过拟合，支撑集里混入了噪声原子'); else disp('残差在噪声范围内，支撑集可信'); end这个判据有效的前提是支撑集大小远小于观测维度 m。如果支撑集本身就接近 m，最小二乘会精确拟合出零残差，这时把残差和噪声能量对比没有意义。把这段检查逻辑加在你的 omp_matlab 主循环之后，比对单纯的信噪比参数要可靠得多，尤其适合批处理仿真里需要自动丢弃失败结果的情况。
本文还有配套的精品资源，点击获取
