简介本资源是一套面向信号处理、机器学习及通信方向本科生与研究生的压缩感知CS算法实践教学包聚焦稀疏信号重建核心问题助力理解奈奎斯特采样之外的高效重构范式。压缩包共7个MATLAB源文件.m涵盖STOMP、SWOMP、SP、IHT、GOMP、OMP和BP七种主流重建算法代码结构清晰、注释完整分别实现随机阈值选择、加权原子筛选、L1范数优化、硬阈值迭代等关键机制便于对比不同策略在噪声鲁棒性、收敛速度与重构精度上的差异。资源体积仅12KB轻量易部署适合课程实验、课程设计及科研入门阶段快速验证理论。目前已有101人学习下载读者可直接运行各算法灵活调整测量矩阵、稀疏度与信噪比参数观察残差演化与支撑集识别过程深入掌握压缩感知从建模到求解的完整技术链条。 前阵子从同事那里拿到一个压缩感知算法实现包名字就叫“压缩感知算法实现.rar”。解压进去看了一眼代码量不算大但真正跑通、跑出稳定结果比我想象中花了更多时间。折腾完以后我发现这份实现其实很有代表性它把稀疏表示、测量矩阵设计、重构算法三块内容完整串了起来而且踩坑点非常典型。这篇文章就记录一下我把它从解压到调优跑通的全过程也聊聊压缩感知里那些代码文档不会写清楚的细节。这份实现适合谁看如果你正在学压缩感知或者做信号处理、图像重建手里有现成数据但不知道怎么套算法那这篇文章能帮你把理论映射到代码。我会讲清每一段代码在干什么、关键参数怎么定、为什么这样定以及我当时是怎么一步步排查问题的。即使你只熟悉Python基础把文中的代码复制下来也可以直接在笔记本上跑一维信号的完整实验。1. 压缩感知项目到底在解决什么问题1.1 一个RAR包背后的核心思路很多人第一次接触压缩感知会把注意力放到“压缩”两个字上。但这份代码包真正的核心其实是在说一件事能不能用比传统奈奎斯特采样少得多的测量次数恢复出原始信号。传统采样讲究采样率必须达到信号最高频率的两倍以上否则混叠。压缩感知换了个思路如果信号在某个变换域里是稀疏的也就是大部分系数为零或接近零那么我可以用一个随机测量矩阵把高维信号投影到低维空间得到一组远小于信号长度的测量值。之后再通过求解一个稀疏优化问题从这组“看起来不完整”的测量里恢复出原始信号。这份RAR里的代码做的就是后面这件事。它没有依赖现成的压缩感知库而是从底层实现了核心算法。我当时看到代码的第一反应是稀疏表示、测量矩阵、重构算法这三样东西是否都被正确处理了。因为只要有一步搞错重构结果就会变成一堆噪声。1.2 拆开看三个关键模块稀疏基、测量矩阵、重构器整个压缩感知流程可以拆成三个模块。第一个是稀疏基。信号本身可能不是稀疏的但它在某个正交基下可以是稀疏的。例如自然图像在小波基下、音频信号在DCT基或短时傅里叶基下大部分系数都接近零。代码里通常用DCT矩阵或傅里叶矩阵来充当这个稀疏基。第二个是测量矩阵。它负责把原始信号从N维压到M维其中M远小于N。这个矩阵有自己的讲究必须满足RIP性质或者具备低相干性。最简单的办法就是生成一个随机高斯矩阵或随机伯努利矩阵这种矩阵在大概率上性质都不错。第三个是重构算法。拿到M维测量值之后怎么恢复出N维信号OMP、基追踪、ISTA、CoSaMP都是常见方法。这份代码包里的重点是OMP它简单直观而且对中等规模问题效果很稳定。这三个模块互相影响我在实际调参时体会很深。稀疏基没有选好重构质量立刻下降测量矩阵归一化没做OMP迭代会出现奇奇怪怪的数值问题重构算法的停止条件设置不对结果要么过拟合噪声要么恢复不完整。1.3 这个项目适合谁参考如果你正在做传感器数据压缩、图像压缩采样、雷达信号处理或者只是想把压缩感知的论文公式变成能跑的代码这份实现都值得研究。它没有复杂的工程依赖核心算法就几十行非常适合用来建立“理论到代码”的映射。我当时的目标很明确用这份代码做一维稀疏信号恢复并演示测量数M对恢复质量的影响。后面我会把完整实验流程和参数计算过程写出来照着做就能复现。2. 压缩感知算法实现的核心模块解析2.1 稀疏表示为什么信号不能直接测如果信号本身不稀疏压缩感知的优势就发挥不出来。这一点很多人会忽视。代码里构造的输入信号如果不是人工构造的稀疏信号那必须先用稀疏基做变换再在稀疏域里进行测量。举个例子一段长度为1024的音频信号它在时域里看起来完全不稀疏。但经过DCT变换之后可能只有前几十个系数是有明显幅度的其余都接近零。这时候我们可以把DCT系数看成稀疏向量alpha原始信号x Psi * alpha。测量过程实际测得的是y Phi * x Phi * Psi * alpha其中Phi是测量矩阵Psi是稀疏基矩阵。这份代码里最需要注意的地方就是感知矩阵A Phi * Psi。很多人在OMP实现里直接拿Phi当作感知矩阵用结果重构系数时一团糟。只有当信号本身已经在稀疏域里表示时才能把Phi直接当感知矩阵。如果不是就必须先计算A Phi Psi再拿A去参与OMP迭代。DCT矩阵的构造方式也不复杂。对于一维信号长度为NDCT基的每个元素可以写成特定余弦函数的值。实际工程中我习惯直接用scipy.fft.dct生成变换矩阵或者用numpy手动构造正交DCT矩阵。手动构造时要小心正交归一化否则重构再逆变换回去时幅度会漂移。2.2 测量矩阵设计M取多少才算“少测”测量数是整个压缩感知实验里最关键的参数。M太小信息量不足M太大压缩感知的意义就没了。理论上有RIP条件但工程上更常用的是经验公式M C * K * log(N / K)。其中K是信号在稀疏域里的非零系数个数N是信号长度C是一个常数通常取2到4。我当时做一个N1024、K30的实验算出log(1024/30)大约等于3.53。按C3来算M最少要大于106。为了留足余量我选了M120大约只有信号长度的12%。这个比例对于一维稀疏信号来说已经能恢复得不错。如果把K调大到50M的下限会变成大约165这时候用M200比较安全。下面是我常用的一组参数对照方便大家快速选初始值信号长度N稀疏度K理论最小MC3建议初始M102410约71100102430约106150102450约166200409650约1992564096100约265350测量矩阵本身用什么分布高斯随机矩阵每个元素独立同分布服从均值为0、方差为1/M的正态分布。伯努利随机矩阵每个元素以独立等概率取±1/sqrt(M)。从效果上看这两种矩阵在中小规模问题上差别不大。高斯矩阵更通用伯努利矩阵在某些硬件实现里更方便因为乘法变成了加减法。RIP条件听起来很数学本质上就是要求测量矩阵不会把两个不同稀疏信号映射到同一组测量值。随机矩阵以很高概率满足这个性质所以代码里直接用随机数生成不需要手动构造确定性矩阵。但有一点要注意每次跑实验固定随机种子否则你会得到不同的随机矩阵实验结果也会不稳定。2.3 重构算法选型OMP、基追踪、ISTA怎么选重构算法是这份代码的核心。常见的四类算法各有优缺点算法思路优势劣势适用场景OMP每次找相关性最大的列逐步扩充支撑集最小二乘更新实现简单速度快需要已知或估计稀疏度K中等规模信号K相对较小基追踪BP求解L1范数最小化理论保证强对噪声更稳通常需要调用优化库速度慢信号较稀疏且量级不大ISTA迭代软阈值收缩实现简单适合大规模收敛慢参数需要调大规模问题GPU上有优势CoSaMP每次选择多列再修剪支撑集比OMP更稳健代码复杂度高一些支撑集变化较快的情形这份实现里为什么会选OMP一方面它足够短便于理解每一步是在干什么另一方面OMP对一维稀疏信号的重构效果已经很能说明问题。如果后面要处理含噪数据再考虑基追踪或者加正则项的ISTA会更好。OMP每次迭代做的事我简单归纳一下先计算感知矩阵A的所有列与当前残差的内积找出绝对值最大的那一列把它加入支撑集然后通过最小二乘求解当前支撑集下的系数再用新的系数计算残差。重复这个过程直到迭代次数达到稀疏度K或残差小到可以忽略。这里有个容易被忽略的细节A的列要归一化。如果测量矩阵做了归一化A的列模长就是1那么内积的大小可以直接比较。如果列模长相差很大OMP选择支撑集时会被大模长的列带偏导致恢复结果不对。这种坑我在后面问题排查部分详细说。3. 实操过程把压缩感知代码跑通3.1 解压后的代码结构怎么看拿到RAR压缩包后我先把文件全部解压出来用tree命令看了眼目录结构。典型实现一般包含这样几个文件核心算法文件omp.py、bp.py稀疏基构造文件basis.py测量矩阵构造文件sensing.py主演示脚本demo.py实验图像或信号数据文件我当时看到这个包的代码第一件事不是急着运行而是逐个文件浏览一遍确认测量矩阵、稀疏基、感知矩阵之间的关系。这个步骤很重要因为很多实现会把稀疏基藏在某个函数里如果你直接调用重构函数不清楚它内部帮你做了哪一步后面出问题很难定位。建议你也按这个顺序读代码先看主脚本怎么调用再看重构函数输入输出最后再看数据生成部分。这样即使代码有命名混乱也能抓住主干。3.2 从零写一个OMP重构函数这一小节给出我调试后最常用的OMP实现。代码不长但每一个细节都关系到成败import numpy as np from numpy.linalg import lstsq def omp(A, y, K, tol1e-6): OMP算法用于稀疏信号恢复 A: M x N 感知矩阵 y: M 维测量向量 K: 稀疏度预估非零系数个数 tol: 残差相对阈值 返回: N 维稀疏系数向量 M, N A.shape x np.zeros(N, dtypecomplex) r y.copy() idxs [] A_conj A.conj().T # 预先计算共轭转置提升循环效率 for _ in range(K): proj A_conj r pos int(np.argmax(np.abs(proj))) if pos in idxs: # 防止重复选择同一列通常意味着残差已经无法提供新信息 break idxs.append(pos) A_sub A[:, idxs] # 最小二乘求当前支撑集下的系数 x_sub, _, _, _ lstsq(A_sub, y, rcondNone) r y - A_sub x_sub if np.linalg.norm(r) tol * np.linalg.norm(y): break x[idxs] x_sub return x那个if pos in idxs的判断很多人会漏掉。理论上理想情况下不会重复选同一列但数值计算里残差可能已经很小相关性最大的仍然是之前选过的列这时候如果继续迭代最小二乘会变得不稳定。我实际调试时这个判断避免了好几次无穷大或NaN的问题。另外A_conj A.conj().T是一种优化。如果A是实数矩阵conj().T就是普通的转置。但如果信号或者稀疏基是复数的比如用到傅里叶基这里就必须取共轭转置否则内积计算结果不对。这一点在第四节还会重点讲。3.3 一维信号完整实验流程有了OMP函数接下来就可以构造完整实验。我用的流程是先生成一个长度为N的稀疏信号alpha在K个随机位置放置非零值然后构造DCT稀疏基Psi构造高斯测量矩阵Phi计算感知矩阵A Phi Psi测量y Phi Psi alpha最后调用omp恢复alpha再通过x Psi alpha恢复原始信号。关键步骤代码np.random.seed(42) N 1024 K 30 M 150 alpha np.zeros(N, dtypenp.float64) support np.random.choice(N, K, replaceFalse) alpha[support] np.random.randn(K) * 1.0 # DCT稀疏基正交形式 Psi np.zeros((N, N), dtypenp.float64) for k in range(N): for n in range(N): if k 0: Psi[n, k] np.sqrt(1.0 / N) else: Psi[n, k] np.sqrt(2.0 / N) * np.cos(np.pi * (2 * n 1) * k / (2.0 * N)) # 高斯测量矩阵 Phi np.random.randn(M, N) / np.sqrt(M) A Phi Psi y Phi x_original # 等价于 A alpha alpha_hat omp(A, y, K) x_hat Psi alpha_hat print(恢复误差, np.linalg.norm(x_hat - x_original) / np.linalg.norm(x_original))运行之后我测试的相对误差通常在1e-8以下。这时候就很清楚地看到只要能正确估计稀疏度KOMP恢复得相当干净。这里有个实验技巧如果你想观察M对恢复质量的影响可以写一个循环固定N和K让M从60逐步增加到200然后记录相对误差。运行完会看到明显的“相位转变”现象M小于某个阈值时误差极大超过阈值后误差迅速下降。这个阈值位置对照前面那个经验公式你会发现基本吻合。弄明白这个曲线你对压缩感知的信息本质会有更直观的理解。3.4 二维图像的分块压缩感知扩展一维信号跑通后我把它扩展到了二维图像。图像信号直接在全局范围内做压缩感知矩阵规模太大N256x256就有65536维直接构造测量矩阵不现实。常规做法是分块处理把图像切成8x8或16x16的小块每一块独立做压缩感知。我当时用8x8块做实验每个块展平成长度64的向量稀疏度K设置为10左右每块测量数M设置为24压缩比接近37.5%。恢复时逐块OMP重构再拼回去。这个方案跑起来很快一块也就毫秒级。分块会有一个副作用就是块与块之间可能出现明显边界。缓解方法有三种一是用重叠块相邻块重叠几个像素恢复后加权平均二是用更大的块比如16x16减少块数但同时增加稀疏度三是恢复后做一次去块滤波。实际项目中重叠块的效果最稳。4. 实测中踩过的坑与问题排查4.1 稀疏度K设多少、重构效果差怎么办OMP最大的参数瓶颈就是K。K设小了支撑集不够恢复结果丢了成分K设大了算法会把噪声和数值误差也当成有效分量恢复结果出现伪峰。我的建议是先别拍脑袋设K。可以运行一次OMP不设固定迭代次数而是设一个较大的上限然后观察残差随迭代次数的变化。正常情况下残差会在前K次迭代中快速下降之后下降速度变得很平缓。这个转折点就是当前数据有效稀疏度的估计。另一个方法是用“稀疏度扫描”依次尝试K1, 2, 4, 8, 16, 32画一条恢复误差随K变化的曲线选择误差最小的K。这个方法稍微费时但非常直观。实测下来K过小误差大K略大误差还能接受K过大误差开始上升曲线呈现一个明显的低谷。4.2 重构结果全是毛刺感知矩阵和归一化的问题我第一次跑自己的数据时恢复结果全是毛刺完全没法看。排查了一圈发现两个问题。第一个问题感知矩阵用错。我把Phi直接传给了OMP忘了乘稀疏基Psi。如果是人为构造的时域稀疏信号碰巧P等于单位阵那直接传Phi没问题。但如果是图像或音频信号在时域不稀疏必须先用Psi做变换。我当时传错之后OMP恢复出来的是稀疏域的系数再逆变换当然对不上。第二个问题测量矩阵没有归一化。很多代码生成Phi时直接np.random.randn(M, N)没有除以sqrt(M)。这样A的每一列模长不同OMP计算内积时模长大的列天然容易胜出支撑集选择就会偏向某些列恢复结果严重失真。解决办法很简单生成Phi后统一除以sqrt(M)或者显式对A的每一列做归一化。4.3 复数信号或复数变换域里的共轭陷阱如果你处理的是频域数据、通信信号或者雷达回波数据往往是复数。这时候OMP里的内积必须做共轭也就是A.conj().T r不是A.T r。我之前在这个地方摔过一次。实数情况下A.T r和A.conj().T r结果一样所以代码在实数信号上跑得正常。一旦换成复数傅里叶基不取共轭会导致相关性计算错误选出的支撑集是完全错乱的。恢复出来的信号实部和虚部都会有问题。用Matlab的话同样要留意A是共轭转置A.是普通转置OMP中必须用前者。4.4 代码跑得慢几个优化小技巧OMP的时间瓶颈在最小二乘和矩阵乘。N跑到数万维时每次迭代lstsq的代价都不小。我的优化思路有三个。第一预计算A的共轭转置避免循环内重复计算。第二最小二乘可以先用QR分解把A_sub的QR分解在迭代中增量更新而不是每次从头解。第三如果M和K都不大直接用np.linalg.pinv(A_sub)也会比lstsq快一点但数值稳定性差一些谨慎使用。对于二维图像的分块处理还有一个语法层面的优化把分块循环写成向量化或并行。用Python写纯for循环跑几百个块会有些慢可以用concurrent.futures并行处理或者用numba加速。我自己的经验是分块数量多时并行收益明显代码改动也不大。5. 从这份代码包里还能扩展出什么5.1 加噪声之后的稳健性真实数据永远带噪声。测量过程可以建模为y Ax n其中n是高斯白噪声。OMP对噪声比较敏感尤其当稀疏度K估计偏大时会把噪声分量也恢复出来。面对含噪数据我有两个调整建议。一是把停止条件从“固定K次”改成“残差小于阈值”让算法自己决定迭代何时停止。阈值可以设为噪声标准差的若干倍。二是改用基追踪或带L2正则的L1最小化也就是LASSO形式这需要引入一个正则参数lambda。lambda调大了恢复结果偏稀疏但可能欠拟合调小了恢复结果含噪。一般用交叉验证或者L曲线来选。5.2 换一种稀疏基试试不同数据的稀疏基选择差别很大。音频信号用DCT或短时傅里叶图像信号除了DCT小波基通常效果更好雷达信号用傅里叶或Gabor基更合适。如果你发现OMP恢复质量上不去不一定是算法问题很可能是稀疏基选得不好。可以做一个快速比较分别用DCT、小波、傅里叶基对同一组数据进行稀疏化计算前K个最大系数的能量占比。能量占比越高说明这个基越合适。对自然图像来说小波基的稀疏性通常优于DCT所以后续如果扩展到图像优先考虑小波稀疏基。5.3 从一维到实时系统的思考做实时网格化或实时信号处理时压缩感知通常不会拿来做整个系统的完整重建而是和前端采集联合设计。测量矩阵在硬件里可以简化成随机下采样位置后端用OMP或ISTA做在线恢复。分块策略在其中很关键。分块越小单次计算越快但块效应越明显分块越大单次OMP迭代越慢但恢复精度更高。这个平衡需要根据具体硬件耗时来调。我个人的经验是先用8x8或16x16块跑通整个链路再根据实际帧率优化。最后再分享一个小技巧。这份代码里的OMP算法其实可以顺手改成一个简单的压缩感知模拟工具你只需要改稀疏基和信号生成部分就能用来验证不同测量矩阵的恢复效果。我后来很多实验都是在这个基础上加噪声、加不同稀疏度很快就看清了算法边界。对于想深入理解压缩感知的人来说这种“把代码拆开改一改再跑一遍”的方式比看十篇论文都管用。本文还有配套的精品资源点击获取 压缩感知算法实践：OMP重构、稀疏表示与测量矩阵调优全解析
简介本资源是一套面向信号处理、机器学习及通信方向本科生与研究生的压缩感知CS算法实践教学包聚焦稀疏信号重建核心问题助力理解奈奎斯特采样之外的高效重构范式。压缩包共7个MATLAB源文件.m涵盖STOMP、SWOMP、SP、IHT、GOMP、OMP和BP七种主流重建算法代码结构清晰、注释完整分别实现随机阈值选择、加权原子筛选、L1范数优化、硬阈值迭代等关键机制便于对比不同策略在噪声鲁棒性、收敛速度与重构精度上的差异。资源体积仅12KB轻量易部署适合课程实验、课程设计及科研入门阶段快速验证理论。目前已有101人学习下载读者可直接运行各算法灵活调整测量矩阵、稀疏度与信噪比参数观察残差演化与支撑集识别过程深入掌握压缩感知从建模到求解的完整技术链条。 前阵子从同事那里拿到一个压缩感知算法实现包名字就叫“压缩感知算法实现.rar”。解压进去看了一眼代码量不算大但真正跑通、跑出稳定结果比我想象中花了更多时间。折腾完以后我发现这份实现其实很有代表性它把稀疏表示、测量矩阵设计、重构算法三块内容完整串了起来而且踩坑点非常典型。这篇文章就记录一下我把它从解压到调优跑通的全过程也聊聊压缩感知里那些代码文档不会写清楚的细节。这份实现适合谁看如果你正在学压缩感知或者做信号处理、图像重建手里有现成数据但不知道怎么套算法那这篇文章能帮你把理论映射到代码。我会讲清每一段代码在干什么、关键参数怎么定、为什么这样定以及我当时是怎么一步步排查问题的。即使你只熟悉Python基础把文中的代码复制下来也可以直接在笔记本上跑一维信号的完整实验。1. 压缩感知项目到底在解决什么问题1.1 一个RAR包背后的核心思路很多人第一次接触压缩感知会把注意力放到“压缩”两个字上。但这份代码包真正的核心其实是在说一件事能不能用比传统奈奎斯特采样少得多的测量次数恢复出原始信号。传统采样讲究采样率必须达到信号最高频率的两倍以上否则混叠。压缩感知换了个思路如果信号在某个变换域里是稀疏的也就是大部分系数为零或接近零那么我可以用一个随机测量矩阵把高维信号投影到低维空间得到一组远小于信号长度的测量值。之后再通过求解一个稀疏优化问题从这组“看起来不完整”的测量里恢复出原始信号。这份RAR里的代码做的就是后面这件事。它没有依赖现成的压缩感知库而是从底层实现了核心算法。我当时看到代码的第一反应是稀疏表示、测量矩阵、重构算法这三样东西是否都被正确处理了。因为只要有一步搞错重构结果就会变成一堆噪声。1.2 拆开看三个关键模块稀疏基、测量矩阵、重构器整个压缩感知流程可以拆成三个模块。第一个是稀疏基。信号本身可能不是稀疏的但它在某个正交基下可以是稀疏的。例如自然图像在小波基下、音频信号在DCT基或短时傅里叶基下大部分系数都接近零。代码里通常用DCT矩阵或傅里叶矩阵来充当这个稀疏基。第二个是测量矩阵。它负责把原始信号从N维压到M维其中M远小于N。这个矩阵有自己的讲究必须满足RIP性质或者具备低相干性。最简单的办法就是生成一个随机高斯矩阵或随机伯努利矩阵这种矩阵在大概率上性质都不错。第三个是重构算法。拿到M维测量值之后怎么恢复出N维信号OMP、基追踪、ISTA、CoSaMP都是常见方法。这份代码包里的重点是OMP它简单直观而且对中等规模问题效果很稳定。这三个模块互相影响我在实际调参时体会很深。稀疏基没有选好重构质量立刻下降测量矩阵归一化没做OMP迭代会出现奇奇怪怪的数值问题重构算法的停止条件设置不对结果要么过拟合噪声要么恢复不完整。1.3 这个项目适合谁参考如果你正在做传感器数据压缩、图像压缩采样、雷达信号处理或者只是想把压缩感知的论文公式变成能跑的代码这份实现都值得研究。它没有复杂的工程依赖核心算法就几十行非常适合用来建立“理论到代码”的映射。我当时的目标很明确用这份代码做一维稀疏信号恢复并演示测量数M对恢复质量的影响。后面我会把完整实验流程和参数计算过程写出来照着做就能复现。2. 压缩感知算法实现的核心模块解析2.1 稀疏表示为什么信号不能直接测如果信号本身不稀疏压缩感知的优势就发挥不出来。这一点很多人会忽视。代码里构造的输入信号如果不是人工构造的稀疏信号那必须先用稀疏基做变换再在稀疏域里进行测量。举个例子一段长度为1024的音频信号它在时域里看起来完全不稀疏。但经过DCT变换之后可能只有前几十个系数是有明显幅度的其余都接近零。这时候我们可以把DCT系数看成稀疏向量alpha原始信号x Psi * alpha。测量过程实际测得的是y Phi * x Phi * Psi * alpha其中Phi是测量矩阵Psi是稀疏基矩阵。这份代码里最需要注意的地方就是感知矩阵A Phi * Psi。很多人在OMP实现里直接拿Phi当作感知矩阵用结果重构系数时一团糟。只有当信号本身已经在稀疏域里表示时才能把Phi直接当感知矩阵。如果不是就必须先计算A Phi Psi再拿A去参与OMP迭代。DCT矩阵的构造方式也不复杂。对于一维信号长度为NDCT基的每个元素可以写成特定余弦函数的值。实际工程中我习惯直接用scipy.fft.dct生成变换矩阵或者用numpy手动构造正交DCT矩阵。手动构造时要小心正交归一化否则重构再逆变换回去时幅度会漂移。2.2 测量矩阵设计M取多少才算“少测”测量数是整个压缩感知实验里最关键的参数。M太小信息量不足M太大压缩感知的意义就没了。理论上有RIP条件但工程上更常用的是经验公式M C * K * log(N / K)。其中K是信号在稀疏域里的非零系数个数N是信号长度C是一个常数通常取2到4。我当时做一个N1024、K30的实验算出log(1024/30)大约等于3.53。按C3来算M最少要大于106。为了留足余量我选了M120大约只有信号长度的12%。这个比例对于一维稀疏信号来说已经能恢复得不错。如果把K调大到50M的下限会变成大约165这时候用M200比较安全。下面是我常用的一组参数对照方便大家快速选初始值信号长度N稀疏度K理论最小MC3建议初始M102410约71100102430约106150102450约166200409650约1992564096100约265350测量矩阵本身用什么分布高斯随机矩阵每个元素独立同分布服从均值为0、方差为1/M的正态分布。伯努利随机矩阵每个元素以独立等概率取±1/sqrt(M)。从效果上看这两种矩阵在中小规模问题上差别不大。高斯矩阵更通用伯努利矩阵在某些硬件实现里更方便因为乘法变成了加减法。RIP条件听起来很数学本质上就是要求测量矩阵不会把两个不同稀疏信号映射到同一组测量值。随机矩阵以很高概率满足这个性质所以代码里直接用随机数生成不需要手动构造确定性矩阵。但有一点要注意每次跑实验固定随机种子否则你会得到不同的随机矩阵实验结果也会不稳定。2.3 重构算法选型OMP、基追踪、ISTA怎么选重构算法是这份代码的核心。常见的四类算法各有优缺点算法思路优势劣势适用场景OMP每次找相关性最大的列逐步扩充支撑集最小二乘更新实现简单速度快需要已知或估计稀疏度K中等规模信号K相对较小基追踪BP求解L1范数最小化理论保证强对噪声更稳通常需要调用优化库速度慢信号较稀疏且量级不大ISTA迭代软阈值收缩实现简单适合大规模收敛慢参数需要调大规模问题GPU上有优势CoSaMP每次选择多列再修剪支撑集比OMP更稳健代码复杂度高一些支撑集变化较快的情形这份实现里为什么会选OMP一方面它足够短便于理解每一步是在干什么另一方面OMP对一维稀疏信号的重构效果已经很能说明问题。如果后面要处理含噪数据再考虑基追踪或者加正则项的ISTA会更好。OMP每次迭代做的事我简单归纳一下先计算感知矩阵A的所有列与当前残差的内积找出绝对值最大的那一列把它加入支撑集然后通过最小二乘求解当前支撑集下的系数再用新的系数计算残差。重复这个过程直到迭代次数达到稀疏度K或残差小到可以忽略。这里有个容易被忽略的细节A的列要归一化。如果测量矩阵做了归一化A的列模长就是1那么内积的大小可以直接比较。如果列模长相差很大OMP选择支撑集时会被大模长的列带偏导致恢复结果不对。这种坑我在后面问题排查部分详细说。3. 实操过程把压缩感知代码跑通3.1 解压后的代码结构怎么看拿到RAR压缩包后我先把文件全部解压出来用tree命令看了眼目录结构。典型实现一般包含这样几个文件核心算法文件omp.py、bp.py稀疏基构造文件basis.py测量矩阵构造文件sensing.py主演示脚本demo.py实验图像或信号数据文件我当时看到这个包的代码第一件事不是急着运行而是逐个文件浏览一遍确认测量矩阵、稀疏基、感知矩阵之间的关系。这个步骤很重要因为很多实现会把稀疏基藏在某个函数里如果你直接调用重构函数不清楚它内部帮你做了哪一步后面出问题很难定位。建议你也按这个顺序读代码先看主脚本怎么调用再看重构函数输入输出最后再看数据生成部分。这样即使代码有命名混乱也能抓住主干。3.2 从零写一个OMP重构函数这一小节给出我调试后最常用的OMP实现。代码不长但每一个细节都关系到成败import numpy as np from numpy.linalg import lstsq def omp(A, y, K, tol1e-6): OMP算法用于稀疏信号恢复 A: M x N 感知矩阵 y: M 维测量向量 K: 稀疏度预估非零系数个数 tol: 残差相对阈值 返回: N 维稀疏系数向量 M, N A.shape x np.zeros(N, dtypecomplex) r y.copy() idxs [] A_conj A.conj().T # 预先计算共轭转置提升循环效率 for _ in range(K): proj A_conj r pos int(np.argmax(np.abs(proj))) if pos in idxs: # 防止重复选择同一列通常意味着残差已经无法提供新信息 break idxs.append(pos) A_sub A[:, idxs] # 最小二乘求当前支撑集下的系数 x_sub, _, _, _ lstsq(A_sub, y, rcondNone) r y - A_sub x_sub if np.linalg.norm(r) tol * np.linalg.norm(y): break x[idxs] x_sub return x那个if pos in idxs的判断很多人会漏掉。理论上理想情况下不会重复选同一列但数值计算里残差可能已经很小相关性最大的仍然是之前选过的列这时候如果继续迭代最小二乘会变得不稳定。我实际调试时这个判断避免了好几次无穷大或NaN的问题。另外A_conj A.conj().T是一种优化。如果A是实数矩阵conj().T就是普通的转置。但如果信号或者稀疏基是复数的比如用到傅里叶基这里就必须取共轭转置否则内积计算结果不对。这一点在第四节还会重点讲。3.3 一维信号完整实验流程有了OMP函数接下来就可以构造完整实验。我用的流程是先生成一个长度为N的稀疏信号alpha在K个随机位置放置非零值然后构造DCT稀疏基Psi构造高斯测量矩阵Phi计算感知矩阵A Phi Psi测量y Phi Psi alpha最后调用omp恢复alpha再通过x Psi alpha恢复原始信号。关键步骤代码np.random.seed(42) N 1024 K 30 M 150 alpha np.zeros(N, dtypenp.float64) support np.random.choice(N, K, replaceFalse) alpha[support] np.random.randn(K) * 1.0 # DCT稀疏基正交形式 Psi np.zeros((N, N), dtypenp.float64) for k in range(N): for n in range(N): if k 0: Psi[n, k] np.sqrt(1.0 / N) else: Psi[n, k] np.sqrt(2.0 / N) * np.cos(np.pi * (2 * n 1) * k / (2.0 * N)) # 高斯测量矩阵 Phi np.random.randn(M, N) / np.sqrt(M) A Phi Psi y Phi x_original # 等价于 A alpha alpha_hat omp(A, y, K) x_hat Psi alpha_hat print(恢复误差, np.linalg.norm(x_hat - x_original) / np.linalg.norm(x_original))运行之后我测试的相对误差通常在1e-8以下。这时候就很清楚地看到只要能正确估计稀疏度KOMP恢复得相当干净。这里有个实验技巧如果你想观察M对恢复质量的影响可以写一个循环固定N和K让M从60逐步增加到200然后记录相对误差。运行完会看到明显的“相位转变”现象M小于某个阈值时误差极大超过阈值后误差迅速下降。这个阈值位置对照前面那个经验公式你会发现基本吻合。弄明白这个曲线你对压缩感知的信息本质会有更直观的理解。3.4 二维图像的分块压缩感知扩展一维信号跑通后我把它扩展到了二维图像。图像信号直接在全局范围内做压缩感知矩阵规模太大N256x256就有65536维直接构造测量矩阵不现实。常规做法是分块处理把图像切成8x8或16x16的小块每一块独立做压缩感知。我当时用8x8块做实验每个块展平成长度64的向量稀疏度K设置为10左右每块测量数M设置为24压缩比接近37.5%。恢复时逐块OMP重构再拼回去。这个方案跑起来很快一块也就毫秒级。分块会有一个副作用就是块与块之间可能出现明显边界。缓解方法有三种一是用重叠块相邻块重叠几个像素恢复后加权平均二是用更大的块比如16x16减少块数但同时增加稀疏度三是恢复后做一次去块滤波。实际项目中重叠块的效果最稳。4. 实测中踩过的坑与问题排查4.1 稀疏度K设多少、重构效果差怎么办OMP最大的参数瓶颈就是K。K设小了支撑集不够恢复结果丢了成分K设大了算法会把噪声和数值误差也当成有效分量恢复结果出现伪峰。我的建议是先别拍脑袋设K。可以运行一次OMP不设固定迭代次数而是设一个较大的上限然后观察残差随迭代次数的变化。正常情况下残差会在前K次迭代中快速下降之后下降速度变得很平缓。这个转折点就是当前数据有效稀疏度的估计。另一个方法是用“稀疏度扫描”依次尝试K1, 2, 4, 8, 16, 32画一条恢复误差随K变化的曲线选择误差最小的K。这个方法稍微费时但非常直观。实测下来K过小误差大K略大误差还能接受K过大误差开始上升曲线呈现一个明显的低谷。4.2 重构结果全是毛刺感知矩阵和归一化的问题我第一次跑自己的数据时恢复结果全是毛刺完全没法看。排查了一圈发现两个问题。第一个问题感知矩阵用错。我把Phi直接传给了OMP忘了乘稀疏基Psi。如果是人为构造的时域稀疏信号碰巧P等于单位阵那直接传Phi没问题。但如果是图像或音频信号在时域不稀疏必须先用Psi做变换。我当时传错之后OMP恢复出来的是稀疏域的系数再逆变换当然对不上。第二个问题测量矩阵没有归一化。很多代码生成Phi时直接np.random.randn(M, N)没有除以sqrt(M)。这样A的每一列模长不同OMP计算内积时模长大的列天然容易胜出支撑集选择就会偏向某些列恢复结果严重失真。解决办法很简单生成Phi后统一除以sqrt(M)或者显式对A的每一列做归一化。4.3 复数信号或复数变换域里的共轭陷阱如果你处理的是频域数据、通信信号或者雷达回波数据往往是复数。这时候OMP里的内积必须做共轭也就是A.conj().T r不是A.T r。我之前在这个地方摔过一次。实数情况下A.T r和A.conj().T r结果一样所以代码在实数信号上跑得正常。一旦换成复数傅里叶基不取共轭会导致相关性计算错误选出的支撑集是完全错乱的。恢复出来的信号实部和虚部都会有问题。用Matlab的话同样要留意A是共轭转置A.是普通转置OMP中必须用前者。4.4 代码跑得慢几个优化小技巧OMP的时间瓶颈在最小二乘和矩阵乘。N跑到数万维时每次迭代lstsq的代价都不小。我的优化思路有三个。第一预计算A的共轭转置避免循环内重复计算。第二最小二乘可以先用QR分解把A_sub的QR分解在迭代中增量更新而不是每次从头解。第三如果M和K都不大直接用np.linalg.pinv(A_sub)也会比lstsq快一点但数值稳定性差一些谨慎使用。对于二维图像的分块处理还有一个语法层面的优化把分块循环写成向量化或并行。用Python写纯for循环跑几百个块会有些慢可以用concurrent.futures并行处理或者用numba加速。我自己的经验是分块数量多时并行收益明显代码改动也不大。5. 从这份代码包里还能扩展出什么5.1 加噪声之后的稳健性真实数据永远带噪声。测量过程可以建模为y Ax n其中n是高斯白噪声。OMP对噪声比较敏感尤其当稀疏度K估计偏大时会把噪声分量也恢复出来。面对含噪数据我有两个调整建议。一是把停止条件从“固定K次”改成“残差小于阈值”让算法自己决定迭代何时停止。阈值可以设为噪声标准差的若干倍。二是改用基追踪或带L2正则的L1最小化也就是LASSO形式这需要引入一个正则参数lambda。lambda调大了恢复结果偏稀疏但可能欠拟合调小了恢复结果含噪。一般用交叉验证或者L曲线来选。5.2 换一种稀疏基试试不同数据的稀疏基选择差别很大。音频信号用DCT或短时傅里叶图像信号除了DCT小波基通常效果更好雷达信号用傅里叶或Gabor基更合适。如果你发现OMP恢复质量上不去不一定是算法问题很可能是稀疏基选得不好。可以做一个快速比较分别用DCT、小波、傅里叶基对同一组数据进行稀疏化计算前K个最大系数的能量占比。能量占比越高说明这个基越合适。对自然图像来说小波基的稀疏性通常优于DCT所以后续如果扩展到图像优先考虑小波稀疏基。5.3 从一维到实时系统的思考做实时网格化或实时信号处理时压缩感知通常不会拿来做整个系统的完整重建而是和前端采集联合设计。测量矩阵在硬件里可以简化成随机下采样位置后端用OMP或ISTA做在线恢复。分块策略在其中很关键。分块越小单次计算越快但块效应越明显分块越大单次OMP迭代越慢但恢复精度更高。这个平衡需要根据具体硬件耗时来调。我个人的经验是先用8x8或16x16块跑通整个链路再根据实际帧率优化。最后再分享一个小技巧。这份代码里的OMP算法其实可以顺手改成一个简单的压缩感知模拟工具你只需要改稀疏基和信号生成部分就能用来验证不同测量矩阵的恢复效果。我后来很多实验都是在这个基础上加噪声、加不同稀疏度很快就看清了算法边界。对于想深入理解压缩感知的人来说这种“把代码拆开改一改再跑一遍”的方式比看十篇论文都管用。本文还有配套的精品资源点击获取