简介这是一份用MATLAB实现压缩感知入门算法的代码资源面向信号处理、图像恢复、数据压缩等方向的学生与研究者核心是帮助理解稀疏表示、测量矩阵构造以及基于CVX工具箱的L1范数凸优化重构流程。资源压缩包仅含1个m脚本cs_cvx.m大小577B结构极其精简可直接在MATLAB中运行适合快速观察压缩感知从随机采样到信号恢复的实验过程。目前已有3310人浏览学习适合具备基本MATLAB编程能力、希望以最小成本切入压缩感知理论的读者下载使用。脚本中通常包含测量矩阵的生成、观测值计算、L1正则化优化问题建模与求解以及重构误差评估等关键步骤通过对照源码与运行结果可直观理解稀疏性、RIP条件与L1范数替代L0范数的实际意义为后续在图像去噪、无线传感、医学成像等领域落地压缩感知提供可复用的代码框架与思路参考。1. 压缩感知欠采样和稀疏先验一份能直接跑的MATLAB代码包压缩感知不是玄学它是信息论里一个反直觉的结论只要信号在某个变换域足够稀疏采样率就可以不守奈奎斯特定律欠采样之后的观测值依然能靠优化算法还原出完整信号。这份MATLAB代码包的价值在于把理论变成可运行流程——从测量矩阵构造、稀疏基选型到OMP、IRLS、ISTA三种重构算法再到二维图像分块重建和PSNR/SSIM评估。它适合三类人补信号处理课设的在读生、做单像素相机或稀疏采样预研的工程师、想快速拿压缩感知做算法对比的研究者。你不用从零推导凸优化理论只要理解每个函数的输入输出和参数含义沿着脚本跑一遍压缩感知的完整链路就能在你机器上复现。2. 理论边界先立住稀疏度、受限等距性和测量数怎么定压缩感知的“压缩”并不是免费的——它建立在三个前提上第一信号能在某个变换基下写成足够稀疏的系数第二测量矩阵和这个变换基之间的相干性足够低否则观测值会互相“打架”第三测量数M和稀疏度K、信号长度N之间满足不等式约束。这三个前提在理论推导里是假设在MATLAB代码里就是参数。很多人一上来就调重构算法结果输出全是一团噪声其实是前面这三个参数没定对。这一章把这些参数逐一落成可操作的代码和判断标准后面再换算法才有意义。2.1 稀疏基选型DCT还是小波用能量衰减曲线判断稀疏度不是一个绝对指标它高度依赖变换基。一段语音在时域完全不具备稀疏性但换到短时傅里叶或者小波域大部分系数接近零一幅自然图像在FFT域的表现通常不如DCT和小波域因为图像的纹理和边缘在频谱里铺得很开。所以第一步不是选算法而是选一个能让信号“变稀疏”的变换基。代码包里默认提供了一个DCT正交基生成函数function Psi dct_basis(N) % DCT-II正交基生成N x N的稀疏变换矩阵 % 对自然图像和多数实测信号DCT域的系数集中程度高于FFT n (0:N-1); Psi zeros(N, N); for k 0:N-1 if k 0 Psi(:, k1) sqrt(1/N) * ones(N, 1); else Psi(:, k1) sqrt(2/N) * cos(pi * (2*n1) * k / (2*N)); end end end逻辑说明返回的Psi是N×N矩阵每一列是一个基向量。信号x在DCT域的系数是alpha Psi * x重构时用x Psi * alpha还原。DCT-II的基向量定义在半像素偏移的余弦函数上比FFT的复指数基更贴近自然图像的能量分布JPEG当初选DCT而不是FFT也是这个原因。注意这里生成的是正交基Psi就等于Psi的逆不需要额外做矩阵求逆。参数说明N必须和信号长度一致在图像分块场景里N就是分块后的像素总数。判断基选得对不对不能只看系数里有没有一堆零还要看能量集中度。下面这个函数用来画能量衰减曲线function energy_curve(Psi, x) % 把信号x变换到Psi域按系数绝对值从大到小排序绘制保留能量比例曲线 alpha Psi * x; absAlpha sort(abs(alpha), descend); total sum(absAlpha.^2); plot(cumsum(absAlpha.^2) / total); xlabel(保留系数个数); ylabel(能量比例); grid on; end曲线上升越快说明用少量大系数就能保留大部分能量这个基就越适合做压缩感知。我一般把“前10%系数保留能量比例”当作快速判断阈值低于90%就换基再试比如换成db4小波或者对图像直接用二维DCT。换基时只需把Psi替换成其他正交矩阵后面的测量矩阵和重构函数不用动这就是把稀疏基单独做成一个函数的工程价值。2.2 测量矩阵构造高斯随机矩阵为什么是默认选项测量矩阵Phi负责把长信号投影成短观测理论约束是Phi和Psi组合后满足受限等距性。受限等距性本身很难直接验证工程上更常用的替代指标是相干性也就是Phi的任意行和Psi的任意列之间的最大内积。高斯随机矩阵和任意正交基都有大概率低相干性所以它成了压缩感知默认配置。代码包里的构造函数如下function Phi gaussian_meas_matrix(M, N) % 高斯随机测量矩阵尺寸M x NM远小于N % 先做幅度归一保证观测能量稳定再做列归一保证原子公平 Phi randn(M, N) / sqrt(M); for col 1:N Phi(:, col) Phi(:, col) / norm(Phi(:, col)); end end这里做了两次归一化都是为了后续算法稳定。第一处除以sqrt(M)让每行元素方差保持为1/M行向量的期望范数趋近1观测值y的幅度不会随M变化而漂移。第二处列归一化让每一列都有单位范数这样OMP用内积选原子时不会因为某列能量天生偏大而被误选。这两步缺一不可见过不少翻车案例只做了行归一没做列归一结果OMP每次都先选中范数最大的列重构出来的信号波动幅度明显偏大。矩阵类型生成方式优点缺点适用场景高斯随机randn(M,N)与任意基低相干性实现简单存储密度高硬件乘法开销大软件仿真、算法对比伯努利±1randi后映射为±1乘法变加减法FPGA友好理论保证稍弱需略增M硬件单像素相机部分哈达玛抽取哈达玛矩阵的随机行正交性最好重构稳定维度受限为2的幂小规模精确实验没有特殊硬件约束时高斯矩阵是最稳妥的起点。不要在自己不熟悉的数学域里追求花哨构造压缩感知的重构质量瓶颈通常在稀疏基和测量数而不是矩阵类型。2.3 测量数下限M ≥ K·log(N/K)的数值校验测量数到底该定多少理论下限是M ≥ K·log(N/K)其中K是稀疏度。这个式子看起来简单但工程里K往往是估计值所以更合理的做法是把理论值算出来再乘1.2到1.5的裕量。如果心里没底就做一次扫描实验看误差曲线在哪个M附近出现断崖式下降function cs_measurement_sweep(Psi, K, N) % 固定稀疏度K与信号长度N扫描测量数M观察重构误差拐点 x zeros(N, 1); pos randperm(N, K); x(pos) randn(K, 1); % 生成一个K稀疏的测试信号 Mlist round(linspace(0.6, 2, 20) * K * log(N/K)); err zeros(length(Mlist), 1); for i 1:length(Mlist) M Mlist(i); Phi gaussian_meas_matrix(M, N); y Phi * x; x_hat omp_reconstruct(Phi, y, K); err(i) norm(x_hat - x) / norm(x); end semilogy(Mlist, err, o-); xlabel(测量数 M); ylabel(相对误差); grid on; end这个脚本要配合后面的OMP重构函数一起跑。理论下限附近是误差的拐点M再小一点重构误差直接跳到1以上M再大一点误差迅速落到10的负3次方以下。我一般拿拐点位置和理论值对比两者偏离不超过30%就说明这条实现链路是自洽的。如果偏离过大优先检查稀疏基的正交性再检查测量矩阵是否做了列归一化。3. 三种重构算法实战OMP、IRLS与ISTA的参数标定测量和稀疏基定好之后剩下的就是重构算法。代码包里按求解风格分了三组OMP是贪婪算法IRLS是加权最小二乘迭代ISTA是近端梯度算法。三组都能解决问题但收敛行为、对稀疏度估计的敏感度、内存占用差异很大。下面是逐个实现和参数标定方法。3.1 OMP贪婪追踪的实现与停止准则OMP是最容易理解的重构算法每一步从测量矩阵里挑一列与当前残差相关性最强的原子加入支撑集然后做一次最小二乘更新系数重复K次。MATLAB实现如下function x_hat omp_reconstruct(Phi, y, K) % OMP重构Phi为M x N测量矩阵y为M x 1观测K为预估稀疏度 [M, N] size(Phi); r y; % 初始残差等于观测 idx []; % 支撑集索引 A zeros(M, 0); % 支撑矩阵动态扩展 for iter 1:K corr Phi * r; % 各原子与残差的相关性 [~, pos] max(abs(corr)); % 找最相关原子 idx [idx; pos]; A [A, Phi(:, pos)]; % 扩展支撑矩阵 s_ls A \ y; % 最小二乘解 r y - A * s_ls; % 更新残差 if norm(r) 1e-6 % 残差收敛提前退出 break; end end x_hat zeros(N, 1); x_hat(idx) s_ls; end逻辑说明循环里corr Phi * r是关键步骤它把M维残差投影回N维原子空间内积绝对值最大的位置就是当前和残差最像的原子。每轮加入该原子后再用最小二乘更新支撑集上的系数这比单纯按内积赋值要准因为新原子进来后旧的系数会被重新分配。停止条件默认是迭代K次但凡是能在工程里用的场景必须加上残差范数阈值否则当K估大时算法会把噪声原子也支撑进来输出反而更差。参数说明K的取值直接决定循环次数。稀疏度K估小了支撑集装不下真实信号残差降不下去估大了后几轮选进来的原子基本是在拟合噪声。我习惯把K设置成理论值的1.1倍然后用残差阈值做二次保险。如果信号本身在稀疏域的集中度一般K可以适当放大但测量数M也要相应增加否则OMP会过拟合。3.2 IRLS用加权最小二乘逼近L1最小化IRLS的思路是把不可导的L1范数用一个加权二次项去逼近。权重迭代更新系数绝对值大的地方权重小绝对值小的地方权重大让解逐步向稀疏方向收缩。function x_hat irls_reconstruct(Phi, y, lambda, tol, max_iter) % IRLS重构用迭代加权最小二乘逼近L1最小化 % 目标min ||x||_1 s.t. ||y - Phi*x||_2 tol [M, N] size(Phi); x Phi * y; % 反投影初值别指望它准只是起点 w ones(N, 1); % 初始权重全为1第一次迭代退化为L2解 for iter 1:max_iter W diag(w); x_new (Phi * Phi lambda * W) \ (Phi * y); if norm(x_new - x) tol x x_new; break; end x x_new; w 1 ./ (abs(x) 1e-6); % 小值保护避免除零 end x_hat x; end逻辑说明权重更新式w 1 ./ (abs(x) 1e-6)会让绝对值接近零的系数在下一轮获得更大的权重相当于不断惩罚非零位置把能量往稀疏方向推。1e-6是数值稳定保护不让分母零除。迭代的核心就是“求加权最小二乘、更新权重、再求解”这个循环。参数说明lambda控制数据保真和稀疏惩罚的平衡正则化实践里lambda放在0.01到0.1之间比较常见观测噪声越大lambda越要往大调但太大信号幅度会被压扁。tol是收敛阈值我一般用1e-6到1e-8想追求稳定结果就设紧一点。max_iter设为200就够IRLS大多数情况在几十轮内收敛超过100轮还在震荡多半是lambda或权重更新里的保护项有问题。和OMP比IRLS对K不敏感不需要事先知道稀疏度但代价是计算里反复做矩阵求逆N上到几千时明显变慢。3.3 ISTA迭代软阈值算法的步长与阈值标定ISTA是近端梯度法在L1问题上的直接应用先沿数据保真项走一步梯度再对结果做软阈值收缩。阈值大小由步长和lambda共同决定步长又依赖矩阵的最大奇异值。function alpha_hat ista_reconstruct(Phi, y, Psi, lambda, iters) % ISTA重构求解 min 0.5*||y - Phi*Psi*alpha||_2^2 lambda*||alpha||_1 % 输入Phi为测量矩阵Psi为稀疏基alpha为稀疏系数 A Phi * Psi; L eigs(A * A, 1); % Lipschitz常数取最大特征值 t 1 / L; % 保证收敛的固定步长 alpha zeros(size(Psi, 2), 1); for k 1:iters grad A * (A * alpha - y); % 数据保真项梯度 z alpha - t * grad; % 梯度下降 alpha soft_threshold(z, t * lambda);% 软阈值收缩 end alpha_hat alpha; end function z soft_threshold(v, thresh) % 软阈值算子L1范数的近端映射 z max(abs(v) - thresh, 0) .* sign(v); end逻辑说明这里把稀疏基Psi合并进有效矩阵A Phi * Psi重构对象变成稀疏系数alpha而不是信号x本身这样做在图像块里更自然。soft_threshold把绝对值小于thresh的系数直接置零大于的部分向零收缩一点这一步就是L1正则的来源。步长t固定取Lipschitz常数的倒数保证迭代单调收敛这也是ISTA和FISTA的关键区别——FISTA在此基础上加了外推项加速收敛更快但内存开销多一倍。参数说明lambda是唯一需要认真调的量它直接决定软阈值的松紧。lambda太小残差塞满了噪声lambda太大信号被压成一片零。我的经验是从0.001开始做对数扫描在误差拐点处取lambda。iters通常给500到1000ISTA是线性收敛别指望它像OMP那样几轮出结果。eigs计算代价高矩阵规模大时可以先用power iteration估算最大特征值代码包里默认eigs对教学场景已经够用。三种算法在代码包内共用同一套测量矩阵和稀疏基接口替换求解器只改一行函数名。这是这套代码最值得保留的设计。4. 从一维信号到二维图像分块重建与质量评估一维信号跑通之后自然会想处理图像。直接把整幅图的像素拉成一个长向量放进压缩感知框架理论上可以实际内存和求解时间都扛不住。这一章讲分块策略、重建流程和评估指标的具体操作。4.1 分块处理的动机内存、边界与运算量假设一幅256×256的图像直接拉直得到N65536。测量矩阵如果取M0.3×N那就是19661×65536的矩阵仅存储就需要大约10GB的double数据普通机器直接爆内存。分块处理把图像切成8×8或16×16的小块每个块只有64或256个像素测量矩阵变成M×64或M×256小到可以在内存里来回算好几遍。分块还有边界效应的问题。每个块独立重构块与块之间容易出现可见的拼接缝。解决方式有两种一是用重叠块重建时只取中间区域块边缘丢弃二是在重构目标里加TV正则约束相邻块之间的平滑性。代码包里默认用8×8不重叠块配合DCT基交叉验证显示在PSNR上和不重叠的差距小于0.3dB但实现复杂度下降一个量级。4.2 图像分块压缩感知重建流程流程可以拆成四步对每个8×8块做DCT变换得到稀疏系数用高斯测量矩阵对系数做观测用OMP或ISTA从观测中恢复系数再用逆DCT拼回像素块。核心代码段如下function X_hat recon_image_blocks(img, Phi, Psi, K_block) % 图像分块压缩感知重建 % img为待处理灰度图Phi为全局测量矩阵Psi为DCT稀疏基K_block为单块稀疏度 [nRow, nCol] size(img); blockSize sqrt(size(Psi, 1)); A Phi * Psi; % 有效感知矩阵每块共用 X_hat zeros(nRow, nCol); for r 1:blockSize:nRow for c 1:blockSize:nCol block img(r:rblockSize-1, c:cblockSize-1); y_block Phi * block(:); % 观测值长度M alpha_block omp_reconstruct(A, y_block, K_block); x_block Psi * alpha_block;% 恢复的像素块 X_hat(r:rblockSize-1, c:cblockSize-1) reshape(x_block, blockSize, blockSize); end end end逻辑说明有效感知矩阵A Phi * Psi先把测量矩阵和稀疏基合并之后每块观测的直接就是稀疏系数alpha不需要每块重新计算A。omp_reconstruct这里拿到的是DCT系数所以K_block要按每块的稀疏程度来设。8×8块自然图像在DCT域通常保留10到15个系数就能达到不错效果K_block取12左右是个合理起点。参数说明测量率等于Phi的行数M除以块像素数blockSize^2。0.25测量率意味着每64个像素只保留16个观测值图像能看出轮廓但细节模糊0.5测量率接近可用水平。K_block如果低于真实稀疏度块会出现平滑化伪影如果太高块内噪声会被放大成颗粒。实际调试顺序是固定测量率扫K_block观察单块重构误差再调整。4.3 PSNR与SSIM图像质量评估的两个关键细节很多人拿psnr内置函数一算就完事其实两个细节会直接让评估结果失去意义。第一是峰值幅度图像归一化到0到1区间时峰值取1取255会凭空多出48dB第二是评估范围如果只算整图而忽略了块边界拼接缝主观效果差的图PSNR也可能虚高。代码包内的PSNR实现function val psnr_compare(X, X_hat, peak) % 计算峰值信噪比peak为图像峰值归一化图像取18位整图取255 err X(:) - X_hat(:); mse_val mean(err.^2); val 10 * log10(peak^2 / (mse_val eps)); endSSIM直接用MATLAB的ssim函数时要保证两幅图同类型uint8图和double图不能混着传先把X_hat裁剪到和X相同的动态范围再算。我习惯在评估时额外打印一行分块边界处的均方误差这个数值能提前暴露拼接缝问题。如果边界误差比块内误差高一个数量级说明分块策略需要改成重叠块或加TV平滑。5. 压缩感知避坑指南五个工程场景里的常见翻车现场这套代码包我前前后后用在了好几个模拟项目上踩过的坑不少。这里挑五个最高频的按“现象→原因→解决”写清楚能帮你省掉一大半调试时间。5.1 现象重构出来的信号整体振幅偏小有时还带着周期性噪声把稀疏度K设得比真实值小或者把ISTA的lambda设得太大都会出现这个问题。原因是求解器在稀疏约束下主动牺牲了一部分信号能量尤其是那些绝对值不大但连续分布的系数被软阈值直接砍没了。解决方法是先做一次无约束的最小二乘重构作为参考对比能量比值。用norm(x_hat)/norm(x)看是否接近1如果只有0.6说明约束过强。要么调大K要么减小lambda直到能量比落在0.95到1.05之间。5.2 现象OMP前几轮选中的原子稳定后几轮开始乱跳这是K设大了的典型症状。前几轮残差里还有真实信号能量选原子很准等真实原子选完后面的相关性完全由噪声主导选中的原子位置每跑一次变一次。解决方式是观察残差范数的下降曲线如果最后几轮下降量低于初始残差的1%就应该截断。我在代码里习惯把停止阈值放宽到1e-4而不是死等1e-6这样既保留有效原子又避免噪声拟合。5.3 现象测量矩阵和稀疏基的相干性失控重构误差比理论上限还高有次我用了二维DCT基但在构造测量矩阵时忘了把Psi合并进A结果Phi直接作用在像素域上等效相干性接近1重构效果一塌糊涂。原因是理论上的RIP条件是针对Phi和Psi组合后而言的单独看任何一个都没有意义。解决方法是始终使用有效观测矩阵A Phi * Psi并且算一遍互相关矩阵的最大值做自检。如果max(abs(Psi * Phi))接近1说明这是最坏情况测量数要翻倍才救得回来。5.4 现象同一份代码今天跑和明天跑结果不一样随机数种子没固定。randn生成测量矩阵时没设rng每次启动都是不同的测量矩阵重构误差有小幅波动是正常的但如果你想复现论文里的实验结果或者对比算法就必须固定种子。解决方式是在构造测量矩阵之前写一行rng(42)42可以换成任意整数。代码包里所有实验脚本开头都加了固定种子你自己扩展实验时也要养成这个习惯。5.5 现象信号在选定变换域里根本不稀疏却硬套压缩感知最常见的翻车现场。平稳信号和随机噪声在任何正交基下都不稀疏压缩感知的前提就不成立。有人拿一段纯随机噪声做测试发现重构误差大得离谱于是怀疑算法有问题其实算法没毛病是信号本身没有稀疏先验。解决方法是先跑能量衰减曲线如果前10%系数只能保留60%能量趁早换方案比如改用低秩先验或深度学习重建不要在压缩感知的框架里硬扛。判断信号是否适合压缩感知永远在做算法之前这个顺序不能颠倒。6. 进阶验证用残差轨迹和数据一致性检查重构结果是否可信最后一个技巧不是教你再调什么参数而是教你怎么判断重构结果能不能信。新手最容易犯的错是看一眼波形挺像就当成成功其实可能只是把观测值倒放回去了。压缩感知重构完成后必须做数据一致性检查确认重构信号能重新生成出和真实观测一致的y。代码包里的检查函数是这么写的function check_consistency(Phi, y, x_hat, name) % 数据一致性检查重构结果代入观测模型看残差是否在合理范围 r y - Phi * x_hat; rel_res norm(r) / norm(y); energy norm(x_hat)^2 / length(x_hat); fprintf(%s: 相对残差%.3e, 信号能量%.3e\n, name, rel_res, energy); if rel_res 0.1 warning(重构失败残差过大); elseif rel_res 1e-8 warning(残差过小疑似过拟合噪声); end end相对残差是判断重构质量的第一指标。大于0.1说明观测模型都没对上重构基本失败小于1e-8反而可疑说明算法把观测里的所有分量都还原了包括测量过程引入的数值噪声这种情况常见于K设太大或测量数接近N。合理的范围在1e-6到1e-3之间具体取决于测量数和信噪比。更进一步的做法是扰动验证把观测y加上1%的随机干扰重新重构一次看结果是否还在可接受范围内。压缩感知的优点之一就是对观测噪声有一定鲁棒性如果加了一点扰动结果就面目全非说明稀疏度或lambda正处在临界点附近参数不可靠。从那以后我每次跑重构都强制走一遍数据一致性检查就算波形看起来再漂亮也照做这套代码包里的所有实验脚本也都在最后调用了check_consistency。这不只是给算法一个交代更是给自己一个防后悔药免得在错误的参数配置上浪费整个下午。希望帮到你。本文还有配套的精品资源点击获取