简介：这份资源面向数字图像处理、计算机视觉与医学成像方向的学习者和研究者，提供二维图像压缩感知重构算法的完整MATLAB实现，帮助读者从代码层面理解CS理论如何用远少于像素数的测量值重建图像。压缩包共26个文件，以18个m脚本为核心，覆盖测量矩阵生成、图像采样、ISTA/FISTA等重构算法与性能评估流程，另含7个fig结果图和1个bmp测试图像，整体约1.59MB，结构紧凑便于逐文件调试。已有4295人学习下载，说明其在教学与实验场景中具备一定参考价值。读者可据此对比高斯、伯努利等不同测量矩阵的重构质量，调整松弛因子与阈值参数，并借助MSE、PSNR指标和可视化结果验证算法效果，也可作为向视频压缩、医学影像、遥感分析及深度学习结合方向拓展的实验起点。

1. 二维图像压缩感知重构：一份能直接跑通的 MATLAB 代码包

做图像处理的人迟早会撞上压缩感知这个坎。采样率明明压到了 30% 甚至更低，重构出来的图却要跟原图肉眼难辨，这中间的落差就是重构算法要填的坑。这份资源给的不是论文里的伪代码，而是一套完整的 MATLAB 程序代码，专门针对二维图像做压缩感知重构。它把测量矩阵生成、稀疏基构造、重构求解、PSNR 对比这几步全串起来了，拿到手改个图片路径就能出结果。适合正在做图像处理大作业、想验证 OMP 或 CoSaMP 这类算法实际效果的人，也适合需要快速搭一个重构 baseline 再往上改的研究场景。MATLAB 版本不用追新，R2016b 之后的都能跑，核心函数没用到太偏的工具箱。

2. 压缩感知重构的数学骨架：从测量矩阵到稀疏基的选型逻辑

2.1 为什么二维图像不能直接套一维压缩感知

一维压缩感知的模型是 y = Φx，x 是长度 N 的列向量，Φ 是 M×N 的测量矩阵，M 远小于 N。但二维图像是 N×N 的矩阵，如果直接按列拉伸成 N²×1 的向量，测量矩阵的尺寸会膨胀到 M×N²，内存和计算量都吃不消。常见做法是分块处理，把图像切成 32×32 或 64×64 的小块，每块单独做压缩感知，最后拼回去。另一种做法是用可分离的测量矩阵，行方向和列方向各乘一个矩阵，把二维运算拆成两次一维运算。这份代码走的是分块路线，块大小在参数区可以改，默认给的是 32×32，兼顾了重构质量和速度。

分块之后有个细节容易被忽略：块与块之间的边界效应。如果每块独立重构再硬拼，块边缘会出现明显的方格感。代码里用了一个简单的重叠分块策略，块间重叠 4 个像素，重构时对重叠区域做加权平均。这个处理不复杂，但有没有它，视觉效果差很多。

2.2 测量矩阵的选取：高斯矩阵还是伯努利矩阵

测量矩阵要满足 RIP 条件，实际中没法直接验证，所以选的都是大概率满足 RIP 的随机矩阵。高斯矩阵是最常用的，每个元素服从 N(0, 1/M)，跟绝大多数稀疏基都不相关。伯努利矩阵元素取 ±1/√M，硬件实现时更友好，但重构质量在高压缩比下会略逊于高斯矩阵。

代码里默认用的是高斯矩阵，生成方式如下：

% 生成 M x N 的高斯测量矩阵，M 为测量数，N 为信号长度 function Phi = generate_gaussian_matrix(M, N) % 每个元素服从 N(0, 1/M)，保证列向量近似正交 Phi = randn(M, N) / sqrt(M); % 对每一列做归一化，避免数值量级差异过大 Phi = Phi ./ repmat(sqrt(sum(Phi.^2, 1)), M, 1); end

这里做了列归一化，原因是后续重构算法对矩阵列向量的能量敏感。如果不归一化，某些列的范数偏大，会导致求解时权重失衡，重构出来的图像局部过亮或过暗。归一化之后每列能量为 1，求解稳定性明显提升。参数 M 就是测量数，M/N 就是采样率，代码里通过采样率反算 M，比如采样率 0.3、块大小 32×32 时，N=1024，M=307。

2.3 稀疏基的选择：DCT 还是小波

自然图像在 DCT 基下是近似稀疏的，在小波基下更稀疏。DCT 基构造简单，就是余弦函数采样，计算快；小波基需要滤波器组，多一层分解，但压缩感知重构质量通常更好。代码里两种基都给了，通过一个 flag 切换。

DCT 基的构造代码：

% 构造 N x N 的 DCT 正交基 function Psi = generate_dct_basis(N) Psi = zeros(N, N); for k = 0:N-1 for n = 0:N-1 if k == 0 Psi(k+1, n+1) = 1 / sqrt(N); else Psi(k+1, n+1) = sqrt(2/N) * cos(pi * k * (2*n+1) / (2*N)); end end end end

这个双重循环在 N=1024 时会比较慢，实际跑的时候可以用 dctmtx 函数替代，MATLAB 自带，一行就能生成。代码里保留手写版本是为了让读者看清 DCT 基的每个元素怎么来的，实际使用时建议换成 dctmtx。小波基用的是 sym4 小波，通过 wavedec2 做两层分解，再把系数拉成列向量。注意小波分解后的系数长度跟原信号长度一致，但排列顺序会影响稀疏性，代码里做了按能量排序的处理，把大系数排在前面，这样重构时阈值收缩更有效。

3. 重构算法实战：OMP、CoSaMP 与 IRLS 的代码实现与参数调优

3.1 OMP 算法的 MATLAB 实现与停止条件设置

OMP 是正交匹配追踪，思路是每次迭代选一个跟残差最相关的原子，然后对已选原子做最小二乘更新。它的优点是实现简单、速度快，缺点是需要知道稀疏度 K，K 估不准就重构不好。

代码里的 OMP 核心循环：

function x_rec = omp_reconstruct(y, Phi, K) % y: 测量向量 M x 1 % Phi: 测量矩阵 M x N % K: 稀疏度 [M, N] = size(Phi); residual = y; support_set = []; x_rec = zeros(N, 1); for iter = 1:K % 计算残差与各原子的内积，找最大相关原子 correlations = abs(Phi' * residual); [~, idx] = max(correlations); support_set = [support_set, idx]; % 在已选原子集上做最小二乘 Phi_sub = Phi(:, support_set); x_sub = pinv(Phi_sub) * y; % 更新残差 residual = y - Phi_sub * x_sub; % 残差能量低于阈值时提前退出 if norm(residual) < 1e-6 break; end end x_rec(support_set) = x_sub; end

停止条件有两个：迭代次数达到 K，或者残差能量低于 1e-6。实际跑的时候，如果 K 设得偏大，OMP 会把噪声也选进来，重构图像出现颗粒感。K 设得偏小，细节丢失。经验值是 K 取 N 的 10% 到 20%，具体看图像纹理复杂度。代码里 K 是通过稀疏度估计函数自动算的，用 DCT 系数排序后取累计能量 95% 的位置作为 K，比固定值靠谱。

3.2 CoSaMP 算法的回溯机制与参数影响

CoSaMP 跟 OMP 的区别在于每次迭代选 2K 个原子，然后回溯剪枝，只保留 K 个。这个回溯机制让它在稀疏度估计不准时更鲁棒，代价是每次迭代计算量翻倍。

代码里 CoSaMP 的关键步骤：

function x_rec = cosamp_reconstruct(y, Phi, K) [M, N] = size(Phi); residual = y; support_set = []; x_rec = zeros(N, 1); max_iter = 50; for iter = 1:max_iter % 步骤1：形成代理信号，找 2K 个最大相关原子 proxy = Phi' * residual; [~, sort_idx] = sort(abs(proxy), 'descend'); new_support = sort_idx(1:min(2*K, N)); % 步骤2：合并支撑集 support_set = union(support_set, new_support); % 步骤3：在合并支撑集上做最小二乘 Phi_sub = Phi(:, support_set); x_sub = pinv(Phi_sub) * y; % 步骤4：剪枝，只保留 K 个最大系数 [~, keep_idx] = sort(abs(x_sub), 'descend'); keep_idx = keep_idx(1:min(K, length(keep_idx))); support_set = support_set(keep_idx); x_sub = x_sub(keep_idx); % 步骤5：更新残差 residual = y - Phi(:, support_set) * x_sub; if norm(residual) < 1e-6 break; end end x_rec(support_set) = x_sub; end

max_iter 默认 50，实际跑的时候大多数图像在 20 次以内就收敛了。如果跑到 50 次还没收敛，通常是测量矩阵条件数太差或者稀疏基选得不合适。回溯剪枝那一步是 CoSaMP 的核心，它保证了支撑集大小始终不超过 K，避免 OMP 那种选错原子就一路错到底的问题。代价是每次迭代要做两次排序，N 大时排序开销不可忽略。

3.3 IRLS 算法：用迭代重加权最小二乘逼近 L0 范数

IRLS 不走贪婪路线，它用一系列加权 L2 问题去逼近 L0 范数。每次迭代根据当前解的幅度调整权重，幅度小的系数权重加大，逼着它们往零走。

function x_rec = irls_reconstruct(y, Phi, max_iter, epsilon) [M, N] = size(Phi); x_rec = Phi' * y; % 初始解用匹配滤波 weights = ones(N, 1); for iter = 1:max_iter % 构造加权矩阵 W = diag(weights); % 求解加权最小二乘：min ||Wx||_2 s.t. Phi*x = y % 用拉格朗日乘子法转为无约束问题 x_rec = W \ (Phi' * ((Phi / W) * Phi' \ y)); % 更新权重，epsilon 防止除零 weights = 1 ./ (abs(x_rec) + epsilon); % 权重变化很小时退出 if norm(weights - 1./(abs(x_rec)+epsilon)) < 1e-4 break; end end end

epsilon 是个关键参数，默认给 1e-3。太大则稀疏性不够，太小则数值不稳定。IRLS 的优点是重构精度通常比 OMP 高，缺点是要解线性方程组，每次迭代计算量比 OMP 大一个量级。图像块小的时候（32×32）还能接受，块大到 128×128 就明显慢了。代码里给 IRLS 设了最大迭代 30 次，实际跑下来 15 次左右权重就基本稳定了。

4. 避坑与排查：重构图像质量不达标时先查这五处

4.1 重构图像全黑或全白

现象：跑完代码，imshow 出来的图要么全黑要么全白，PSNR 低到个位数。

原因：测量矩阵没有归一化，或者稀疏基的逆变换用错了。全黑通常是重构系数全为零，全白是系数爆炸。检查 generate_gaussian_matrix 里有没有做列归一化，以及 DCT 基的转置有没有用对。

解决：在重构完成后加一句 x_rec = x_rec / max(abs(x_rec(:)))，先归一化再显示。如果还是不对，把测量矩阵换成 dctmtx 的前 M 行，这个是确定性的，能排除随机矩阵的问题。

4.2 PSNR 卡在 20dB 上不去

现象：采样率提到 0.5 了，PSNR 还是 20 出头，跟论文里的 35dB 差很远。

原因：稀疏度 K 设得不对，或者分块大小太小。32×32 的块在采样率 0.5 时，M=512，对纹理复杂的块来说测量数还是不够。另外 OMP 的 K 如果按固定值给，遇到平滑块会过拟合，遇到纹理块会欠拟合。

解决：把块大小改成 64×64，采样率不变的情况下每块测量数翻倍。K 改成自适应，用 DCT 系数累计能量 95% 来定。如果还不行，换 CoSaMP，它对 K 的估计误差容忍度更高。

4.3 重构时间随图像尺寸指数增长

现象：512×512 的图跑了十分钟还没出结果，内存占用飙升。

原因：没有分块，直接把整图拉成 262144×1 的向量，测量矩阵尺寸 131072×262144，内存直接爆。或者分块了但每块都重新生成测量矩阵，重复计算。

解决：确认代码走的是分块路径，块大小不超过 64×64。测量矩阵生成一次就存下来，所有块共用。如果内存还是紧张，把测量矩阵改成稀疏存储，MATLAB 的 sparse 函数能省不少内存。

4.4 块边界出现方格状伪影

现象：重构图像放大看，块与块之间有明显的亮度跳变，像马赛克。

原因：分块时没有重叠，每块独立重构后硬拼，块边缘的系数不连续。

解决：启用重叠分块，重叠像素数设为块大小的 1/8 到 1/4。重构后对重叠区域做加权平均，权重用余弦窗，中间高两边低。代码里 overlap 参数默认是 4，块大小 32 时相当于 12.5% 重叠，够用了。

4.5 换图片后重构完全失效

现象：用代码自带的测试图跑得好好的，换成自己的图就重构不出来。

原因：图片尺寸不是块大小的整数倍，分块时边界处理出错。或者图片是彩色三通道，代码只处理了灰度。

解决：读图后先做尺寸裁剪，把长宽都裁成块大小的整数倍。彩色图先用 rgb2gray 转灰度，或者对三个通道分别做压缩感知再合并。代码里加一句 size_check 函数，自动裁剪并给出提示。

5. 从重构到评估：PSNR、SSIM 与采样率曲线的自动化脚本

跑通单张图的重构只是第一步，真正做实验的时候需要批量跑不同采样率、不同算法，然后画对比曲线。这份代码里带了一个评估脚本，把 PSNR 和 SSIM 的计算封装好了，直接调用就行。

PSNR 的计算：

function psnr_val = calc_psnr(original, reconstructed) % 确保数据类型一致 original = double(original); reconstructed = double(reconstructed); mse = mean((original(:) - reconstructed(:)).^2); if mse == 0 psnr_val = Inf; return; end max_val = max(original(:)); psnr_val = 10 * log10(max_val^2 / mse); end

注意 max_val 取的是原图最大值，不是固定的 255。如果图像归一化到 [0,1] 了，max_val 就是 1，PSNR 数值会跟 [0,255] 范围差 20dB 左右，对比论文数据时要统一量纲。

SSIM 的计算用 MATLAB 自带的 ssim 函数，但要注意它默认是彩色图模式，灰度图要加参数：

ssim_val = ssim(reconstructed, original, 'DynamicRange', max(original(:)));

批量实验的脚本结构是三层循环：外层采样率从 0.1 到 0.5 步进 0.1，中层算法遍历 OMP、CoSaMP、IRLS，内层遍历测试图片。每次跑完把 PSNR 和 SSIM 存到矩阵里，最后用 plot 画曲线。

我一般会加一个计时器，记录每种算法的平均重构时间。因为实际选算法的时候不只看质量，还要看能不能在可接受的时间内跑完。OMP 在采样率 0.3 时单块重构大概 0.05 秒，CoSaMP 是 0.12 秒，IRLS 要 0.8 秒。整图 512×512 分 256 个块，OMP 总共 13 秒左右，IRLS 要三分钟以上。如果做实时性要求高的场景，IRLS 基本不用考虑。

还有一个技巧是固定随机种子。测量矩阵是随机生成的，不同种子跑出来的 PSNR 会有 0.5dB 左右的波动。做对比实验时把 rng(42) 写在循环外面，保证所有算法用的是同一个测量矩阵，这样对比才公平。这个细节论文里经常不写，但实际复现时不知道就会觉得结果对不上。

从那以后我每次跑压缩感知实验，都先把随机种子固定、量纲统一、块大小和重叠参数确认一遍，这三步走完再开始跑数据。希望帮到你。

本文还有配套的精品资源，点击获取