简介稀疏表示是信号处理与机器学习中的核心方法通过少量非零元素逼近原信号可显著揭示数据内在结构广泛用于图像去噪、压缩感知MRI重建与字典学习等场景。这份RAR压缩包共6个M文件、体量仅4KB代码量虽少但流程闭环主程序、PCA降维模块、L1范数优化求解器、样本读取、精度计算与预处理等环节衔接清晰覆盖从数据准备、稀疏编码到性能评估的完整实验链路。资源面向机器学习初学者和需要快速复现稀疏表示实验的研究者与作者博文阐述的LASSO、BPDN及K-SVD等经典算法配套可直接在MATLAB环境中运行调试帮助读者从代码层面理解稀疏编码与字典学习的核心实现思路。目前已有6565人学习下载适合信号处理、图像分析方向学习者作为入门例程对照研究也可作为课程作业或论文复现的参考代码。1. 稀疏表示到底在解决什么问题一个让 MATLAB 代码从「能跑」到「会用」的转折点你在 MATLAB 命令行里敲过x A \ b觉得线性方程不过如此。但稀疏表示想做的事完全相反它不追求“唯一解”而是在无穷多个可行解里找出那个系数向量绝大多数元素都等于零的解。这个思路在信号压缩、图像去噪、特征提取里反复出现也是压缩感知的前置基础。很多人第一次接触稀疏表示是被“内含完整 MATLAB 代码”的标题吸引进来结果卡在字典构造和迭代算法上。这篇笔记就把整套路径铺平从数学定义到完整可跑的 MATLAB 工程再到参数怎么调、噪声来了怎么救全程用手写代码说话。适合正在做算法对比或写论文的从业者也适合从「调用工具箱」走向「看懂算法」的 MATLAB 用户。2. 稀疏表示的数学骨架字典、稀疏系数与 OMP 的收束逻辑要落地一套 MATLAB 代码先得把“稀疏表示”这四个字拆成能写进矩阵的东西。它不像深度学习那样有一堆网络层核心数学模型其实只有一行等式但这一行等式背后的求解逻辑决定了你在 MATLAB 里写的是循环还是调用现成函数。2.1 稀疏表示的语言从线性组合到 l0 范数约束假设你有一个信号x维度是N×1。稀疏表示说的是这个信号可以被一个过完备字典D维度N×M且M N里的少量列向量线性组合出来。这里的每一列叫作“原子”M个原子张成一个比你原始信号空间更大的空间所以同一个x会有无数种表达方式。数学上写成x D * alphaalpha是M×1的系数向量。稀疏表示的目标是让alpha里的非零元素尽可能少也就是最小化||alpha||_0l0 范数统计非零元素个数。为什么非得少因为少意味着你抓住了信号的本质结构。比如一段音频在某个字典下可能只需要 5 个原子的线性组合就能重建剩下的几百个系数全是零那这 5 个原子的位置和权重就是这段音频的“指纹”。这个思路用在去噪上尤其直观噪声在字典下通常是分散的需要很多原子才能凑出来而干净信号只需要少数几个原子。所以当你强加“稀疏”这个约束时噪声自然就被丢掉了。2.2 l0 问题的直接解是组合爆炸OMP 为什么成为工程首选理论上最完美的求解方式是直接最小化 l0 范数。但现实很残酷要从M个原子中选K个需要遍历C(M, K)种组合。当M256、K6时组合数已经超过十亿量级MATLAB 根本算不动。这就是为什么实际工程里没人硬解 l0。工程上有两条路一条是凸松弛把 l0 换成 l1 范数Basis Pursuit用线性规划或cvx工具箱求解另一条是贪婪算法最常见的就是正交匹配追踪OMP。OMP 的哲学很朴素每次选一个和当前残差最相关的原子选完就用最小二乘更新系数然后算新的残差重复 K 次。它不保证全局最优但实际效果在字典原子相关性不高时非常接近最优解。我一般首选 OMP原因很实际20 行 MATLAB 就能写完不需要额外工具箱迭代过程每一步都可以打印出来看排错非常方便。相比之下凸松弛方法虽然理论更优美但求解器一介入就变成黑匣子出了问题很难查。2.3 字典是稀疏表示的真正变量三种常见字典的选型理由很多人以为稀疏表示的核心是算法其实算法只是搬运工字典才是决定效果上限的变量。同一段信号在 A 字典下稀疏度可能是 5在 B 字典下稀疏度可能变成 50效果天差地别。常见的字典有三类字典类型构造方式适合场景注意点过完备 DCT 字典用离散余弦变换的基向量在不同频率上采样列数大于行数平滑信号、音频、图像块原子间相关性随过完备率上升K 不能设太大小波字典用小波基函数构造需要 Wavelet Toolbox非平稳信号、边缘丰富的图像工具箱依赖重小波基的选择影响大随机高斯字典生成高斯随机矩阵再列归一化压缩感知理论验证理论性质好但物理意义弱不利于解释对入门项目过完备 DCT 字典是最稳妥的选择代码自己就能写不需要任何工具箱而且 DCT 对大多数自然信号都有不错的稀疏性。等你把 OMP 跑顺了再换成 K-SVD 学习字典也不迟。3. 用 MATLAB 跑通最小可复现的稀疏表示工程完整代码与每一步解释这一章直接给一套能跑的完整工程。我习惯把它拆成三个文件一个主脚本负责生成测试数据和验证结果一个my_omp.m函数实现核心算法。这样做的原因是函数单独放文件里MATLAB 的报错信息会精确到行列排错体验比混在一个脚本里舒服得多。3.1 先构造一块能产生“真稀疏”的测试数据做算法验证最怕的是“假数据”。如果你拿一段随机噪声去跑稀疏表示结果毫无意义。所以我先构造一份“真稀疏”的系数再用字典合成信号这样我们心里有底真实的稀疏度是几、非零位置在哪。% test_sparse.m - 稀疏表示主脚本 % 参数区信号长度 N字典原子数 M真实稀疏度 K_true N 128; % 信号长度 M 256; % 字典原子数过完备字典M N K_true 6; % 真实稀疏度 % 1. 构造过完备 DCT 字典 k (0:N-1); % 行索引N×1 j (0:M-1); % 列索引1×M % 经典 DCT-II 形式M 列对应 M 个离散基函数 D sqrt(2/N) * cos(pi/N * (k 0.5) * j); % 第 0 列直流分量的归一化系数要单独调整 D(:, 1) D(:, 1) / sqrt(2); % 列归一化每个原子的能量归一为 1这一步对 OMP 极其关键 col_norm sqrt(sum(D.^2, 1)); D D ./ col_norm; % 2. 生成一份“真稀疏”的系数 alpha_true zeros(M, 1); alpha_true([13 47 99 132 201 240]) [0.8 -0.6 1.2 0.5 -0.9 0.7]; % 3. 用字典把信号合成出来 x D * alpha_true;逻辑说明k和j构造了一个二维采样网格cos项生成的是离散余弦基。M比N大所以这个字典有 256 个原子而信号只有 128 维属于过完备字典。第 0 列对应直流基DCT-II 的定义里它需要额外乘1/sqrt(2)。列归一化必须做否则 OMP 在计算内积时会偏向能量大的原子造成系统性选错。参数说明N128是信号长度对应你实际业务里的信号维度M256是字典原子总数经验上取2N到4N之间效果较好。alpha_true里那 6 个非零位置和值是“真值”后面恢复出来的结果要和它对比才能检验算法是否正确。注意 MATLAB 里取列可以用D(:, [13 47 99])一次取出多列这在调试稀疏系数时非常好用。3.2 手写 20 行 OMP 函数核心迭代与最小二乘更新OMP 的函数实现是整个工程的心脏。我见过很多版本把while循环写得花里胡哨但核心就四步找最相关原子、更新已选原子集合、最小二乘求解、更新残差。% my_omp.m - 正交匹配追踪 % 输入D - N×M 归一化字典 % x - N×1 观测信号 % K - 最大迭代次数稀疏度上限 % 输出alpha - M×1 稀疏系数 function alpha my_omp(D, x, K) N size(D, 1); M size(D, 2); alpha zeros(M, 1); r x; % 残差初始为信号本身 selected zeros(1, K); % 记录选中的原子索引 A zeros(N, K); % 已选原子组成的子字典 for iter 1:K % 计算残差与所有原子的内积绝对值最大者最相关 corr D * r; [~, idx] max(abs(corr)); % 防止同一个原子被重复选中浮点误差下可能发生 while any(selected(1:iter-1) idx) corr(idx) 0; [~, idx] max(abs(corr)); end selected(iter) idx; A(:, iter) D(:, idx); % 在已选原子张成的空间里做最小二乘 coef A(:, 1:iter) \ x; % 残差 信号 - 已解释部分 r x - A(:, 1:iter) * coef; % 残差已小到机器精度继续迭代没有意义 if norm(r) 1e-12 break end end % 回填稀疏系数只有被选中的位置非零 alpha(selected(1:iter)) coef; end逻辑说明每次迭代先算D * r这是所有原子与当前残差的内积内积绝对值越大说明这个原子和残差方向越一致。选中最相关的原子后把它加入A然后用A(:, 1:iter) \ x求解当前最优系数。这里的核心逻辑是用反斜杠而非inv因为最小二乘解要求解一个可能是病态的系统反斜杠在数值稳定性上比显式求逆好得多。残差更新完进入下一轮。参数说明K是最大迭代次数它就是你假定的稀疏度上限。如果K设得比真实稀疏度小恢复就会缺原子设得过大OMP 会把噪声也拟合进来。代码里加了norm(r) 1e-12的提前终止条件这个阈值在无噪声场景下可以保证迭代自动停在正确的稀疏度上。3.3 主脚本串联从字典到信号再到稀疏系数恢复有了测试数据和 OMP 函数主脚本的最后一段就是把这套流程串起来验证恢复结果。% 接 test_sparse.m 已生成的 D 和 x % 这里假设我们知道真实稀疏度先用真实值 K_true 来验证 alpha_est my_omp(D, x, K_true); % 找出恢复系数中非零的位置 idx_est find(abs(alpha_est) 1e-6); idx_true find(abs(alpha_true) 1e-6); % 重建信号并计算归一化重建误差 x_rec D * alpha_est; recon_err norm(x - x_rec) / norm(x); fprintf(真实稀疏度%d恢复稀疏度%d\n, length(idx_true), length(idx_est)); fprintf(原子位置一致%s\n, mat2str(sort(idx_est))); fprintf(重建误差%.2e\n, recon_err);逻辑说明find(abs(alpha_est) 1e-6)是提取稀疏系数里非零位置的标准写法。由于浮点计算会有微小误差理论上应该为 0 的位置可能残留 1e-15 级别的数值所以用一个阈值卡掉。重建误差norm(x - x_rec) / norm(x)是归一化指标小于 1e-12 说明恢复几乎完美。如果一切正常你会看到恢复的原子位置和真实的[13 47 99 132 201 240]完全一致重建误差在 1e-13 量级。这个结果说明代码逻辑没有错。接下来才能放心地加噪声、调参数、做对比实验。3.4 三个可控实验无噪恢复、有噪恢复、稀疏度不足工程验证要一次准备三个场景无噪声验证算法正确性加噪声验证鲁棒性稀疏度不足验证边界行为。% 实验2有噪声场景 x_noisy x 0.01 * randn(N, 1); alpha_noisy my_omp(D, x_noisy, 8); % 稀疏度上限放宽到 8 x_denoised D * alpha_noisy; % 计算去噪前后的信噪比单位 dB snr_in 20 * log10(norm(x) / norm(x_noisy - x)); snr_out 20 * log10(norm(x) / norm(x_denoised - x)); fprintf(输入SNR%.2f dB重建SNR%.2f dB\n, snr_in, snr_out); % 实验3稀疏度不足 alpha_poor my_omp(D, x, 3); % 只给 3 次迭代 x_poor D * alpha_poor; fprintf(稀疏度不足时重建误差%.2e\n, norm(x - x_poor) / norm(x));逻辑说明实验 2 把稀疏度上限放宽到 8是因为噪声会分散能量给一点余量才能让 OMP 把真实原子找全。实验 3 故意只给 3 次迭代这时重建误差会明显变大因为 6 个真实原子只找回了 3 个。这两个实验能帮你建立直觉K 设大会拟合噪声K 设小会丢失信号结构。参数说明0.01是噪声标准差约等于信号能量的 1% 量级。8这个值是按经验取的等于真实稀疏度加 2给噪声留出缓冲。实际项目中噪声标准差的估计可以通过std(x)或中值估计法得到不要拍脑袋定。4. 稀疏表示落地必调的参数稀疏度 K、字典规格与停止条件的边界同样的代码在不同参数下结果可能从“完美恢复”变成“完全乱套”。这一章把三个关键参数的边界说清楚避免你把 OMP 当成黑盒调参。4.1 稀疏度 K 不是拍脑袋用残差曲线定 K真实场景里你往往不知道信号的稀疏度是多少。这时候最可靠的办法是画残差曲线让 K 从 1 递增到一个较大的值观察重建残差的变化。% 残差曲线分析脚本 K_max 30; history zeros(K_max, 1); for k 1:K_max a my_omp(D, x, k); history(k) norm(x - D * a); end % 画图观察拐点 plot(1:K_max, history, o-); xlabel(迭代次数 K); ylabel(重建残差 ||x - D\alpha||); grid on;逻辑说明当 K 小于真实稀疏度时每增加一个原子残差会大幅下降当 K 超过真实稀疏度后新增原子只能解释微小噪声或冗余成分残差下降明显变缓。这个“肘部拐点”就是信号在该字典下的有效稀疏度。这个手段不需要任何先验知识比拍脑袋定 K 靠谱得多。参数说明K_max取 30 还是 50取决于你对信号复杂度的估计。如果残差曲线一直平滑下降没有明显拐点说明这个字典选得不好信号在字典下根本不稀疏这时候调 K 没有意义回头换字典才是正路。4.2 字典规格 N×M 怎么定过完备率与原子相关性的边界字典的行数 N 由信号维度决定能动的是列数 M。M 越大字典表达力越强但原子之间的相关性也会上升OMP 会开始“分不清”两个相似的原子。M 取值过完备率原子互相关系数实际效果M N1 倍完备正交基0正交稀疏度可能较高表达力受限M 2N2 倍约 0.1~0.3平衡点入门首选M 4N4 倍约 0.3~0.6表达力强但 OMP 易选错原子M 8N8 倍接近 0.8不推荐OMP 基本失效互相关系数的定义是字典任意两列内积绝对值的最大值即mu max(|D(:,i) * D(:,j)|), i ~ j。可以用一行 MATLAB 算出来G D * D; G abs(G - diag(diag(G))); mu max(G(:));参数说明当mu超过 0.5 时OMP 的行为就开始不稳定因为残差在相似原子上的投影几乎一样。如果你的字典mu偏高可以改用稀疏正则化方法如 Basis Pursuit或者做一步 Gram 矩阵的阈值收缩但入门阶段最简单的是降低过完备率。4.3 停止条件与噪声门限不要把 K 当唯一选项很多刚接触 OMP 的人把 K 当成唯一参数这是误解。K 只是“最大迭代次数”真正的停止条件有三种达到 K 次、残差低于绝对门限、残差相对能量低于比例门限。无噪声场景下norm(r) 1e-12的绝对门限很可靠。但一旦有噪声这个门限永远达不到OMP 就会一直迭代到 K把噪声全部拟合进来。正确的做法是把停止条件改成相对门限% 带噪声场景下推荐的停止条件 noise_sigma 0.01; % 噪声标准差需要估计 stop_threshold noise_sigma * sqrt(N) * 2; % 残差应降到的水平 relative_threshold 0.05; % 或按信号能量的 5% 作为门限 % 在迭代循环里 % if norm(r) stop_threshold || norm(r) / norm(x) relative_threshold % break % end参数说明stop_threshold的理论依据是当残差接近噪声能量时剩下的成分基本都是噪声继续迭代只会拟合噪声。noise_sigma可以通过diff(x)的中值绝对偏差估计得到这在 MATLAB 里用mad(diff(x))一句话就能算。相对门限0.05是工程经验值适合大多数信号处理场景。5. 常见问题与避坑MATLAB 跑稀疏表示最容易翻车的 5 个点代码不长但每个跑过 OMP 的人都交过学费。这一章把我见过的高频坑按“现象 → 原因 → 解决”写清楚照着排查能省半天时间。5.1 现象脚本里直接定义 function 报错新手最容易踩的坑是把my_omp函数直接写在主脚本中间或末尾然后运行时报错Function definitions are not supported in this context. Functions can only be defined in a script if the script is a function file.原因MATLAB 对脚本文件中函数的位置有严格限制。R2016b 之后脚本末尾可以放局部函数但前提是脚本里不能有普通的脚本代码执行到函数定义处更规范的做法是函数单独存文件。解决把my_omp保存为独立的my_omp.m文件主脚本test_sparse.m放在同一个文件夹下。或把主流程包进一个function test_sparse()里再把my_omp作为局部函数放在文件末尾。我一般用第一种因为独立文件在报错时能直接定位到函数内部某一行。5.2 现象OMP 迭代了几轮后选中的原子重复残差不再下降你打印每次选中的原子索引发现第 4 次选的和第 2 次是同一个。这是因为你构造字典时跳过了列归一化能量大的原子总是被重复选中还有一种可能是浮点误差导致残差在已选原子方向上的分量没有完全清零。原因字典列未归一化或 OMP 代码里没有处理“残差与已选原子内积不为零”的数值误差。解决在字典构造后做D D ./ sqrt(sum(D.^2, 1))确保每列能量为 1。同时在选原子时加一个while去重逻辑把已经选过的索引对应的内积置零后重新选。这两个措施同时做基本不会再翻车。5.3 现象重建误差很小但恢复的稀疏系数和真值完全对不上你把alpha_est的非零位置打印出来发现和alpha_true差了十万八千里但重建误差norm(x - x_rec)却很小你觉得“结果挺好”其实已经错了。原因过完备字典的原子之间存在相关性同一个信号可能有多组稀疏表示。你恢复出来的那一组虽然稀疏但不是真值那一组。这在压缩感知里叫“表示不唯一”。解决区分两个目标——如果你要的是“信号重建”那误差小就够了如果你要的是“恢复真实稀疏系数”必须检查字典的互相关系数mu。当mu大于 0.3 时不要期待系数能和真值完全一致。另一个思路是用 K-SVD 学习字典学习出来的字典原子更匹配信号结构但那是另一个话题。5.4 现象加了噪声后选出来的原子全乱套无噪声时恢复得很漂亮一旦x_noisy里加了0.01 * randnOMP 第一步就选错了原子后面全跟着错。原因OMP 每次选的是和残差“最相关”的原子噪声在高维空间里与某些原子天然有较大内积。当噪声能量不可忽略时第一步就把噪声方向当成了信号方向。解决先估计噪声水平然后把 OMP 的停止条件从“固定 K 次”改成“残差降到噪声水平就停”。更工程化的做法是先对信号做一次轻量平滑比如smoothdata(x_noisy, sgolay, 5)再进 OMP能显著降低噪声对原子选择的干扰。如果噪声很强改用 CoSaMP 或 SP 这类对噪声更鲁棒的贪婪算法。5.5 现象用 inv 和 pinv 求解结果在迭代中被浮点误差放大了有人喜欢写coef inv(A * A) * A * x或者pinv(A) * x。明明数学上等价于A \ x但结果在 OMP 多步迭代后越来越差。原因inv(A*A)需要先算A*A这个矩阵在原子相关性强时会病态条件数很大求逆会放大浮点误差。而反斜杠A \ x用的是 QR 分解或最小二乘求解器数值稳定性好得多。解决统一用coef A(:, 1:iter) \ x。如果矩阵接近秩亏可以用lsqminnorm(A, x, 1e-12)做带容差的最小二乘。记住一个原则MATLAB 里凡是能写成x A \ b的就不要写成x inv(A) * b这不是玄学是数值稳定性。6. 把同一套代码拿去干实事信号去噪验证与 OMP 落地清单前面五章的代码已经能跑通但真正拿到业务里用还需要一套验证方法。我最常做的是把长信号切成短块逐块做稀疏表示去噪然后观察信噪比变化和稀疏度的稳定性。这里有一个快速自检清单——每次用 OMP 做去噪前我都会逐项核对字典是否列归一化mu是否低于 0.4稀疏度 K 是否通过残差曲线选定停止条件是否包含噪声门限求解是否用了反斜杠而非inv。这五个检查点只要有一个不满足结果基本会翻车。我自己的血泪教训是第一次拿 OMP 做音频去噪时狠狠调了一天 K效果始终不如意后来发现是字典列没有归一化OMP 把能量大的原子全选完了真信号反而被丢掉。那次之后我养成了一个习惯所有稀疏表示代码落地的第一件事不是调参而是打印一行常规检查结果fprintf(字典大小%dx%d互相关系数%.3f列能量范围[%.2f, %.2f]\n, ... size(D,1), size(D,2), mu, min(col_norm), max(col_norm));这一行能挡住一大半低级错误。稀疏表示这个方向的价值其实不在那套数学有多深而在你能否用最短的代码把信号拆成“几个原子 一堆零”的样子。把第三章的工程代码跑熟、把第四章的参数边界摸清你就已经超过大多数只会在工具箱里点鼠标的人了。希望帮到你。本文还有配套的精品资源点击获取