简介压缩包提供二维图像压缩感知重构算法的完整MATLAB实现面向信号处理、图像处理方向的研究者、学生及工程师。内容覆盖测量矩阵生成、图像采样、重构算法与性能评估包含OMP、Wavelet_OMP、CS_OMP等典型程序并配有fig可视化文件与示例图像便于对照调试与复现实验。资源共26个文件其中18个m脚本为核心代码7个fig图表用于结果展示另有1个bmp测试图像压缩包大小约1.59MB结构紧凑实用。已有4295人学习下载常用于算法对比、参数调优及课程设计等场景。通过研读代码可掌握迭代软阈值、正交匹配追踪等重构思路并能进一步扩展至不同测量矩阵效果比较、噪声环境下的鲁棒性实验以及医学成像、遥感分析等应用领域是一份兼顾理论理解与动手实践的优质参考资料。1. 二维图像的压缩感知重构一篇能直接跑通的 MATLAB 实战笔记一张 256×256 的灰度图直接按像素采集需要 65536 个测量值压缩感知理论告诉你用其中 5%10% 的随机线性测量配合稀疏重构就能把图像恢复出来。二维图像的压缩感知重构算法程序代码做的就是把这句理论落地成可运行的 MATLAB 工程。这篇笔记从测量矩阵构造、稀疏基选择写起一直到 OMP 迭代、分块策略和参数调法每个关键参数都给出取值依据。适合刚接触压缩感知、想绕开论文公式直接出图的开发者也适合手头代码能跑但重构总翻车、想知道问题出在哪的人。理论只讲必要部分重点是可复现代码和排错思路。2. 压缩感知理论骨架稀疏表示、测量矩阵与分块的工程理由压缩感知重构的起点不是算法而是三个前提信号可稀疏表示、测量矩阵与稀疏基不相关、测量数大于稀疏度若干倍。三个条件都满足欠定方程 yΦx 才有唯一稀疏解。二维图像之所以是压缩感知最经典的应用场景是因为自然图像在 DCT 或小波域天然稀疏这给下面所有代码提供了成立基础。2.1 图像在哪个变换域稀疏DCT 基与稀疏度 K 的估计CS 恢复图像的前提是信号可压缩。对二维图像来说像素域的灰度值天然不平缓直接把像素当 x 去求稀疏解没有意义必须先在某个正交基 Ψ 下做变换x Ψθ其中 θ 只有 K 个显著非零项这个 K 远小于图像维度。这时观测方程变成 y ΦΨθ Aθy 是低维测量A 是 M×N 的传感矩阵。因为 M N方程欠定但若 θ 足够稀疏求解 min‖θ‖₀ s.t. ‖y-Aθ‖≤ε 就能恢复。实际算法用 OMP 这类贪心策略逼近最优解省去组合爆炸的 l0 搜索。图像处理里最常用的稀疏基是 DCT 和小波。DCT 对平滑区域的能量集中度很好16×16 块内几乎全部能量集中在前几十个低频系数小波对边缘和纹理更友好但 MATLAB 里构造可逆小波矩阵要绕一步。工程第一版我一般选 DCTdctmtx 一行生成正交矩阵计算快分块重构时稳定性足够。FFT 基也有应用但周期性延拓在图像块边界会产生振铃不如 DCT 平稳。K 怎么估把一块 16×16 图像做 DCT 后按幅值降序排序计算累计能量前 K 个系数能量占全块能量的比例 ≥95% 就认为 K 有效。多数自然图像在 16×16 块内取 K1530 就能保住 95% 以上能量。如果实测重构结果细节大量丢失优先怀疑 K 取小而不是换基。2.2 测量矩阵为什么选高斯随机RIP 与实现成本测量矩阵 Φ 的作用是把 N 维信号压缩成 M 维测量。理论上要求 Φ 对任意 K 稀疏信号满足约束等距性质对所有 K 稀疏 θ存在常数 δ 使 (1-δ)‖θ‖₂² ≤ ‖ΦΨθ‖₂² ≤ (1δ)‖θ‖₂²保证不同稀疏信号压缩后仍可区分。验证一个矩阵是否满足 RIP 是组合问题工程上没人真去验证直接选高斯随机矩阵。高斯随机矩阵的每个元素独立同分布于标准正态再做 sqrt(M) 归一化。它的优势是任意固定正交基 Ψ 与它组成传感矩阵 A 后以高概率满足 RIP实现只要一行 randn和图像内容不耦合不会因为图像本身有周期性而产生系统性失真。部分傅里叶矩阵虽然利于内存但需要处理索引重排和共轭转置第一版代码不建议碰。我的做法就是固定 rng 后 randn(M, N)/sqrt(M)结果可复现出问题再换。2.3 分块测量全图全局测量为什么会内存爆炸如果对整张 256×256 图像直接构造传感矩阵图像向量长度 N65536取 10% 测量率 M≈6554A 的维度是 6554×65536元素约 4.3 亿个double 存储占 3.4GB这还只是单精度中间量。分块后呢16×16 的块向量长度 256M8 时 A 只有 8×256 个元素几十 KB全图所有块复用同一个 A内存开销可忽略矩阵乘法还能复用。分块代价是块与块之间失去耦合重构后接缝处可能出现方块效应这个坑在第五部分给对策。对首次实现来说分块是把理论跑通的关键决策矩阵小、速度快、内存稳值得先接受块效应再逐步优化。3. 完整 MATLAB 实现分块测量与 OMP 重构的闭环代码3.1 主程序测量矩阵、稀疏基与分块测量% cs_image_omp_demo.m % 二维图像压缩感知重构分块测量 OMP DCT clear; clc; rng(2024); % 固定随机种子保证结果可复现 img imread(cameraman.tif); % 换成你自己的 8 位灰度图 img im2double(img); [rows, cols] size(img); block 16; % 分块大小rows 和 cols 需能被 block 整除 M 8; % 每块测量数决定压缩率 M/block^2 K 24; % 稀疏度估计OMP 迭代上限 D dctmtx(block^2); % block^2 维正交 DCT 矩阵 Psi D; % x Psi * theta Phi randn(M, block^2) / sqrt(M); % 高斯测量矩阵 A Phi * Psi; % 传感矩阵 A Phi * Psi numRow rows / block; numCol cols / block; Y zeros(M, numRow * numCol); % 保存每块的测量向量 idx 1; for r 1:numRow for c 1:numCol xb img((r-1)*block1 : r*block, (c-1)*block1 : c*block); Y(:, idx) Phi * xb(:); % y Phi * x测量端实现 idx idx 1; end end稀疏基用 dctmtx 的转置做反变换矩阵 Psi传感矩阵 A 不需要单独存储全尺寸按块复用。测量循环只做矩阵乘没有隐式维度问题。参数 row/col 通过 numRow/numCol 保证整除不符合时先 imresize 或裁剪。测量向量 Y 的维度是 M × 块数后续重构按同一顺序取列。3.2 OMP 函数原子选择、最小二乘与残差更新function theta omp_solve(A, y, K) % OMP 求 y A * theta e假设 theta 稀疏度为 K [~, N] size(A); At A; theta zeros(N, 1); r y; % 初始残差 selected zeros(K, 1); % 记录被选原子下标 tol_r 1e-6 * norm(y); % 残差相对阈值 for s 1:K corr At * r; % 残差与所有原子的内积 if s 1 corr(selected(1:s-1)) 0; % 禁止重复选同一原子 end [~, pos] max(abs(corr)); % 选内积绝对值最大的原子 selected(s) pos; As A(:, selected(1:s)); % 已选原子构成的子矩阵 theta_ls As \ y; % 最小二乘求当前系数 r y - As * theta_ls; % 更新残差 if norm(r) tol_r % 残差足够小则提前停止 break; end end theta(selected(1:s)) theta_ls; end每次迭代先算残差与所有原子的内积选出相关性最大的原子加入支撑集再用支撑集上的最小二乘获得当前系数最后计算新残差。最小二乘那一步是关键它让已选原子的系数整体重估而不是各自独立求内积因此收敛更快。残差阈值取 norm(y) 的 1e-6 倍避免用绝对阈值在大动态范围图像上误停。K 一旦达到循环停止。OMP 对稀疏度很敏感K 大于真实稀疏度会把噪声也拟合进来K 太小则细节丢。初始 K 设为 M/3 到 M/2 之间M8 时 K24 是偏保守的上限。3.3 重构端分块回溯与质量输出rec zeros(rows, cols); idx 1; for r 1:numRow for c 1:numCol theta omp_solve(A, Y(:, idx), K); xb_rec Psi * theta; % 反变换回像素域 rec((r-1)*block1 : r*block, (c-1)*block1 : c*block) reshape(xb_rec, block, block); idx idx 1; end end rec_disp min(max(rec, 0), 1); % 裁剪到 [0,1] figure; subplot(1,2,1); imshow(img); title(原图); subplot(1,2,2); imshow(rec_disp); title(sprintf(OMP 重构, PSNR %.2f dB, psnr(rec_disp, img)));重构时每块独立调用 omp_solve返回的是 DCT 系数用 Psi 反变换得到像素块再放回原位置。显示前必须裁剪到 [0,1]否则 psnr 会因为越界报错imshow 也会自动截断产生错觉。PSNR 在 2035 dB 之间是正常区间低于 20 dB 先查测量矩阵是否归一化高于 35 dB 基本说明测量率太浪费应该降 M 看实际压缩率。3.4 参数表一组能直接起步的默认值| 参数 | 示例值 | 作用与调整方向 | | block | 16 | 决定向量维度和局部平稳性8 更快32 更依赖图像纹理 | | M | 8 | 每块测量数压缩率 M/block²至少是 K 的 23 倍 | | K | 24 | 稀疏度上限过小丢失细节过大引入噪点 | | rng 种子 | 2024 | 固定测量矩阵便于对比实验 | | tol_r | 1e-6 | 残差相对阈值越小迭代越多 |提示第一版调试不要同时调 M 和 K固定一个参数扫另一个否则很难判断变化来源。4. 算法选型OMP 之外SAMP、CoSaMP 与 TV 各自的适用场景4.1 SAMP稀疏度未知时用步长逼近OMP 最大的弱点是要预先知道 K。实测中 K 对重构质量很敏感K20 和 K30 可能一个模糊一个带噪。SAMP 的思路是以步长 L 分阶段扩展支撑集每阶段先按残差相关性选入 L 个候选原子用最小二乘估计后看残差是否下降如果残差不降反增就进入下一阶段并增大 L。这样稀疏度未知也能收敛代价是迭代次数增加、对步长初值敏感。% SAMP 核心迭代片段简化版 % 输入 A, y, stepSize输出 theta r y; selected []; stage 1; L stepSize; while norm(r) 1e-6 * norm(y) stage * stepSize M [~, pick] sort(abs(A * r), descend); candidate union(selected, pick(1:L)); % 候选集 theta_ls A(:, candidate) \ y; r_new y - A(:, candidate) * theta_ls; if norm(r_new) norm(r) stage stage 1; % 残差没降扩大步长 L stage * stepSize; else selected candidate; r r_new; end end这段是思想验证用的简化版没做剪枝。实际用 SAMP 要注意stepSize 一般取 13太大绕过最优解太小迭代爆炸。适用场景是 K 完全不清楚、需要自动定阶的探索阶段。4.2 CoSaMP用回溯剪枝换稳定性CoSaMP 每次迭代选 2K 个原子进候选集再通过最小二乘和剪枝只保留其中贡献最大的 K 个。它比 OMP 多了一次裁剪抗噪声能力更强而且每轮都重新估计全部已选原子的系数不会像 OMP 那样一旦选错原子就一路错到底。% CoSaMP 单次迭代核心示意 % r 为残差selected 为上一轮支撑集 corr A * r; [~, idxTop] sort(abs(corr), descend); candidate union(selected, idxTop(1:2*K)); % 选 2K 个候选 theta_ls A(:, candidate) \ y; [~, keep] sort(abs(theta_ls), descend); selected candidate(keep(1:K)); % 保留贡献最大的 K 个 r y - A(:, selected) * theta_ls(keep(1:K));CoSaMP 适合稀疏度已知但噪声偏大的场景。对二维图像来说如果 OMP 重构后背景噪点明显而细节还在值得换 CoSaMP 试一轮。它的代价是每轮要解两次最小二乘速度比 OMP 慢约一倍但对初版调试可以接受。4.3 TV 重构面向图像梯度稀疏的先验OMP 家族都在 DCT 域做稀疏但 DCT 对块边界的连续性没有约束。TV全变分模型直接假设图像梯度稀疏目标函数是 min‖x‖_TV s.t. ‖y-Φx‖₂≤ε重构结果天然抑制方块效应和颗粒噪声对医学图像这类大片平滑区域效果好。代价是求解过程是凸优化迭代比 OMP 慢一个数量级而且需要现成的优化求解器不适合只想快速验证算法的场景。我的选型习惯是第一版无脑 OMP如果图像重构后细节够但区域不干净换 CoSaMP如果主要是方块接缝问题优先考虑重叠分块而不是直接上 TV因为 TV 的调参成本更高。4.4 四种算法选型表| 算法 | 稀疏度依赖 | 每轮选原子数 | 抗噪性 | 速度 | 典型场景 | | OMP | 固定 K | 1 | 一般 | 快 | 首次跑通稀疏度有估计 | | SAMP | 自适应 | 步长 L | 中 | 中 | K 未知自动定阶 | | CoSaMP | 固定 K | 2K 剪回 K | 较强 | 较慢 | 噪声明显稳定性优先 | | TV | 不需要 | 无 | 强 | 慢 | 平滑区域多允许离线 |5. 二维图像重构的 5 个常见坑现象、原因与排查5.1 坑一全图构造测量矩阵直接内存崩溃现象代码在构造 A Phi * Psi 或 Phi randn(M, rows*cols) 时卡死任务管理器内存打到上限或 MATLAB 直接报“内存不足”。原因全图向量维数是 rows*cols对 256×256 图像是 65536M 取 10% 时测量矩阵有 4.3 亿个元素double 存储 3.4GB还没算传感矩阵乘积。解决按块测量A 的尺寸降到 M×block²若确实需要全局测量改用部分傅里叶或结构化测量算子不要实例化稠密矩阵。分块是更简单、对算法验证足够的选择先跑通再谈全局。5.2 坑二K 取太大重构图像出现颗粒噪声现象把 K 从 15 调到 40PSNR 不升反降图像上有明显的孤立亮点像撒了盐。原因M8 时测量向量只有 8 维OMP 最多只能可靠恢复约 M/2 个自由度K40 远超过测量信息量算法用额外原子强行拟合测量噪声。解决K 初始取 M/2 上下。我习惯先在小图上扫 K画 PSNR 曲线% 小图快速扫 K避免每轮全图重跑 Klist 4:2:20; pK zeros(size(Klist)); for i 1:numel(Klist) recK cs_reconstruct_64(img64, 16, 8, Klist(i)); pK(i) psnr(recK, img64); end [~, bestIdx] max(pK); K_best Klist(bestIdx);cs_reconstruct_64 是只处理 64×64 小图的封装函数几秒能出结果。K 扫描曲线会在某个拐点后明显回落那个拐点就是当前测量率下的信息量上限。5.3 坑三OMP 重复选原子支撑集条件数恶化现象OMP 循环里连续两次迭代选到同一列最小二乘出现矩阵接近奇异警告残差停滞不降。原因常见于忘记把已选原子从相关投影中排除或者传感矩阵本身两列高度相关内积最大位置重复指向同一列。解决循环里对 selected 下标位置的内积值直接置零代码在 3.2 已经实现。若置零后仍反复选到相近列说明 A 的列相干性太高检查 Phi 是否有 sqrt(M) 归一化Psi 是否方正交矩阵。5.4 坑四分块重构后出现明显方块接缝现象重构图像整体轮廓正常但每个 16×16 块的边缘有一条亮度跳变线均匀区域也呈格子状。原因独立分块重构没有对块边界做一致性约束两块的系数估计在边界处不连续DCT 本身的周期延拓也放大了边界跳变。解决两招。第一招是重叠分块步长取 block/2每个像素被相邻块覆盖重构后重叠区域做加权平均接缝明显软化。第二招是评估时避开边界PSNR 计算时把四周各裁掉 block 像素只统计内部区域避免边界跳变拉低指标。TV 重构本身带平滑先验也能缓解但调参成本高。5.5 坑五PSNR 算出 NaN 或数值离谱现象psnr 返回 NaN或者同一算法两次跑出来一个 18dB 一个 33dB。原因重构值超出 [0,1] 范围psnr 要求 double 输入在 [0,1]不同随机种子导致测量矩阵变化重构质量波动大图像本身不是 [0,1] 范围时 im2double 漏写。解决显示和计算前统一 cliprec_disp min(max(rec, 0), 1); p psnr(rec_disp, img);评估时固定 rng 种子做算法对比才公平。PSNR 波动的正常范围是 ±2dB超过 5dB 先怀疑测量率或 K 参数漂了不是算法坏了。6. 重构质量怎么量化固定种子、PSNR/SSIM 与快速验证流程验证一套实现是否真跑通我习惯按三步走。第一步固定 rng对 64×64 小图做全参数扫描block 固定 16M 从 4 到 16K 取 M/2画 PSNR 曲线曲线应单调上升并趋缓。如果 PSNR 完全不动说明测量或重构链路有逻辑错误不是调参问题。第二步选定一组好的 M/K换三张不同类型的灰度图平滑、纹理、混合看 PSNR 和 SSIM 是否都在合理区间。第三步记录每次实验的 rng 种子和输入参数避免“这次效果好下次复现不了”。封装一个评估函数会更顺手function [p, s] eval_cs(img, block, M, K) % 把前面的主流程封装成函数 cs_reconstruct rec cs_reconstruct(img, block, M, K); rec min(max(rec, 0), 1); p psnr(rec, img); s ssim(rec, img); fprintf(PSNR %.2f dB, SSIM %.3f\n, p, s); endSSIM 比 PSNR 更能反映结构保持程度尤其是判断方块接缝和边缘模糊。SSIM 在 0.9 以上说明重构结构完整0.7 以下就要回查分块策略和稀疏基选择。更细的习惯测量矩阵固定后不要每次都重新生成。把 Phi、Psi、A 作为缓存变量测量率扫描时唯一变的是 MA 需要重新生成但 Psi 可以复用。这样调试时不会因为随机性干扰判断。我踩过的教训是早期没固定 rng同一组参数跑两次 PSNR 差 3dB当时还以为是代码有 bug查了半天后来发现只是测量矩阵变了。固定种子、先小图、再换图这套流程比任何理论推导都更能帮你确认实现对不对。希望帮到你。本文还有配套的精品资源点击获取
