压缩感知稀疏重构算法详解:FOCUSS原理与MATLAB实现
2026/8/31 12:30:02 网站建设 项目流程

简介:本资源是一套面向信号处理与压缩感知领域初学者及科研人员的稀疏重构算法实践代码包,聚焦FOCUSS等主流稀疏求解方法,解决低维测量下高维稀疏信号精确重建这一核心问题,适用于医学成像、雷达信号处理、无线传感等典型CS应用场景。压缩包共含10个MATLAB(.m)源文件,总大小仅9KB,涵盖FOCUSS单/多通道实现(FOCUSS_Single.m、FOCUSS_Multiple.m)、基追踪(BP.m)、正交匹配追踪(OMP_fun.m)、BPDN同伦算法(BPDN_homotopy_function.m)及SBL贝叶斯稀疏学习(SBL_C_fun.m)等关键算法模块,辅以primal/dual/inverse等通用更新函数,结构清晰、模块解耦,便于理解迭代机制与算法差异。目前已有622人学习下载,读者可直接运行验证不同算法在相同仿真条件下的重构精度与收敛特性,快速掌握稀疏建模、残差迭代、L1正则化等核心思想,并为算法改进与工程部署提供可调试的基准实现。

1. 内容整体设计与思路拆解

1.1 这个压缩感知资源包到底解决什么问题

先把这个资源包的老底揭开。标题里挂着"几种常见的稀疏重构算法代码",核心是 FOCUSS,同时也带上了压缩感知和稀疏重构两个大背景词。这类压缩包在网上很常见,但很多人下载回来解压完,打开发现一堆 .m 文件,跑又跑不出来,看又看不懂,最后只能躺在硬盘里吃灰。我写这篇东西,就是想把这类资源包彻底讲透,让拿到代码的你能跑通、能看懂、能改着用。

压缩感知的测量模型就一个式子:y = Ax + n。y 是 M 维观测向量,A 是 M×N 的测量矩阵(M 远小于 N),x 是 N 维原始信号,n 是噪声。如果 x 本身是稀疏的——也就是只有 K 个非零元素,K 远小于 N——那么理论上我们可以从欠定方程 y = Ax 中把 x 恢复出来。这就是压缩感知的核心承诺:采样量远低于奈奎斯特频率,却仍能完美重建。

但"能重建"和"怎么重建"是两回事。从 M 个方程解 N 个未知数,直接求逆是无穷多解的。稀疏重构算法的任务就是在无穷多解中挑出最稀疏的那个。FOCUSS 就是这类算法里比较经典的一条路线,它用迭代加权的方式逐步逼近稀疏解,原理不复杂,代码量也不大,非常适合作为入门压缩感知的第二个算法——第一个通常是 OMP。

这个资源包适合谁?两类人。第一类是刚接触压缩感知的学生,算法原理学过一遍,但不知道代码怎么写,跑起来是什么效果,需要一份能直接跑的参考实现。第二类是做信号处理落地项目的工程师,手里有测量数据,想试试不同稀疏重构算法的效果,需要一个能对比的代码集。不管你是哪类,这篇都能让你把资源包里的东西用起来。

1.2 为什么 FOCUSS 值得单独拎出来讲

OMP 这类贪婪算法有个前提:你要么知道稀疏度 K,要么通过停止条件去猜。但实际工程里 K 往往未知,猜错了结果就很尴尬。FOCUSS 不用预先给定稀疏度,它走的是另一条路——从最小二乘解开始,通过不断调整权重矩阵,把能量逐步集中到少数分量上。

FOCUSS 的数学基础是最小化 lp 范数(0 < p ≤ 1)。为什么是 lp 范数?因为 p=2 的最小二乘解太平滑,p=0 的范数能体现稀疏性但它是非凸的、组合爆炸式的难解,而 0<p≤1 的 lp 范数虽然也是非凸的,但它比 l0 更容易处理,而且能有效逼近稀疏解。FOCUSS 就是通过迭代加权最小二乘(IRLS)的方式,把"最小化 lp 范数"这个目标间接求解出来。简单说,它就是在能量最集中的方向上不断收缩,最终逼出一个稀疏解。

这里要区分一下:BP(Basis Pursuit,基追踪)也做稀疏重构,但它解的是 l1 凸优化,全局最优但计算量大,需要调用线性规划工具箱;OMP 速度快但依赖稀疏度先验;FOCUSS 的定位是介于两者之间——不需要稀疏度,速度比 BP 快,性能在多数情况下接近 l1 方法。这也是为什么很多资源包会把 FOCUSS、OMP、BP 放在一起做对比:三个算法代表了三条完全不同的设计思路,各有各的脾气。

1.3 资源包代码结构的拆解与定位

这类压缩包里出现的文件,按功能大致可以分三组。第一组是核心算法实现,比如 focuss.m、omp.m、bp_magic.m 这类,一个文件对应一个算法,输入是测量矩阵 A、观测 y 和参数,输出是重构出来的稀疏信号 x_hat。第二组是工具函数,比如生成稀疏信号的 gen_sparse.m、计算重构误差的 nmse.m、画对比图的 plot_results.m,这些是辅助你做实验的脚手架。第三组是主脚本,比如 demo_focuss.m、compare_algorithms.m,作用是串起整个流程:造数据、加噪声、跑算法、出指标、画图。

一般拿到压缩包先别急着跑主脚本,我建议按这个顺序来:先打开 focuss.m 看核心实现,理解输入输出;然后检查工具函数是否齐全;最后再跑 demo。实战中很多人一上来就运行 demo,报错了就懵了,其实八成是路径没配对、工具箱缺失、或者 MATLAB 版本不兼容。先把代码结构摸清楚,后面出问题也好定位。

2. 核心细节解析与实操要点

2.1 FOCUSS 算法原理一张纸讲明白

FOCUSS 的迭代式看起来有点吓人,但拆开看就三层:

x_k = W_k · A^T · (A · W_k · A^T)^{-1} · y

其中 W_k = diag(|x_{k-1}|^(1-p/2))。

初始化时,x_0 一般取最小二乘解 A^+ y,也就是 x_0 = A^T (A A^T)^{-1} y。迭代开始后,每一步都用上一次的结果去更新权重矩阵 W,然后把一个加权后的最小二乘问题解出来,得到新的 x_k。权重矩阵的核心逻辑是:上一次解中幅值大的分量,在下一轮迭代中会获得更大的权重,从而进一步放大;幅值小的分量权重被压低,逐步被"挤"向零。这个过程重复下去,解会越来越集中在少数几个非零元素上,直到满足终止条件。

p 值的选择直接影响结果。p 越小,稀疏性越强,但收敛越不稳定;p=0 理论上对应最稀疏解,实际容易震荡;p=1 时稳定性好,近似于 l1 范数优化,但稀疏性略弱。工程上常用 p=0.5 作为折中,我在实际代码里默认给这个值。还有一些改进版本会做正则化,把迭代式变成 x_k = W_k A^T (A W_k A^T + λI)^{-1} y,λ 是正则化参数,这种处理在含噪场景下能显著抑制噪声放大。

2.2 代码包中 FOCUSS 算法的具体实现

这个压缩包里的 focuss.m 核心函数,目前主流的写法基本长这样,我贴出来一份可以对照着看的伪代码:

function x_hat = focuss(A, y, p, max_iter, tol) % A: M*N 测量矩阵 % y: M*1 观测向量 % p: lp norm parameter, e.g., 0.5 % max_iter: 最大迭代次数 % tol: 收敛门限 [M, N] = size(A); % 初始化: 最小二乘解 A_pinv = pinv(A); x = A_pinv * y; for k = 1:max_iter % 更新权重矩阵 w = abs(x).^(1 - p/2); W = diag(w); % 加权最小二乘更新 x_new = W * A' * inv(A * W * A' + 1e-6 * eye(M)) * y; % 收敛判断 if norm(x_new - x) / norm(x) < tol x = x_new; break; end x = x_new; end x_hat = x; end

注意这里加了个 1e-6 的对角微扰,作用是防止 AWA' 奇异导致求逆失败。我第一次实现 FOCUSS 时没加这个,矩阵秩亏的时候直接报错,加了微扰之后稳定性明显提升。稀疏度不依赖外部输入,这是 FOCUSS 相对 OMP 的核心差异。

顺带说一句,这个函数里的 inv() 可以改成 (AWA' + 1e-6*eye(M)) \ y 的形式,求解速度更快,数值稳定性也更好。实际工程中优先用反斜杠运算符,尽量避免显式求逆。

2.3 三个关键参数:p、迭代次数、终止门限

p 值上面说过,我再补充一点:p 的取值范围是 0 到 1 之间,但不要把 p 调得过低(比如 0.1),否则迭代极易发散。我自己测试下来,p=0.5 对于大部分场景是安全的,p=0.8 到 1 适合噪声较大的场景,p 偏小适合无噪或低噪场景。如果你处理的信号稀疏度比较高——比如 10 个非零值分布在 500 维里——p=0.3 到 0.5 能给到更好的稀疏性,但前提是你已经验证了 A 满足某种条件(比如 RIP,受限等距性质)。

迭代次数和经验值有关。我见过有人设 100 次,其实 FOCUSS 在干净场景下通常 20 次左右就收敛了,超过 50 次基本就是 p 太小或者噪声太大。设置 100 次不是不行,但浪费时间。终止门限 tol 我一般设 1e-6,相对变化量低于这个值就认为收敛。在含噪场景下,tol 设太严没有意义,因为噪声本身会阻止解完全稳定,设 1e-4 就够。

还有一个容易被忽略的点:观测 y 的尺度会影响迭代。如果 y 的量纲很大(比如幅值上千),数值计算中矩阵求逆的精度会受影响。建议先对 y 做归一化,迭代完成后再把幅值还原回去。这个细节在多数演示代码里会被省略,但遇到重构结果出现数值异常时,优先检查这个。

3. 实操过程与核心环节实现

3.1 环境准备与数据构造

实操部分我用 MATLAB 来跑,因为压缩感知领域的经典代码基本都是 MATLAB 写的,这个压缩包里的 .m 文件也默认你装了 MATLAB。如果你没有 MATLAB,用 GNU Octave 也能运行大部分代码,只是个别绘图函数需要小改。

第一步是构造一个能测试算法的场景。我们要生成一个 256 维的稀疏信号,稀疏度 K=10,也就是 10 个非零值,其余全是 0。观测维度 M 取 64,测量矩阵 A 用随机高斯矩阵——每个元素独立同分布,均值为 0,方差为 1/M。高斯随机矩阵是压缩感知里最常用的测量矩阵,因为它以高概率满足 RIP 条件。

rng(42); % 固定随机种子,保证实验可重复 N = 256; K = 10; M = 64; x_true = zeros(N, 1); % 随机挑 K 个位置,赋随机幅值 idx = randperm(N, K); x_true(idx) = randn(K, 1) * 5; % 高斯随机测量矩阵 A = randn(M, N) / sqrt(M); % 无噪声观测 y = A * x_true;

这里 randn(K,1)5 让非零元素幅值大于 1,方便后面观察 FOCUSS 的权重集中效果。测量矩阵除以 sqrt(M) 是为了归一化列的能量,这是一种常见做法,让 A 各列的期望范数保持稳定。如果不做这个归一化,AA' 的特征值分布受 M 影响较大,重构性能会有波动。

3.2 运行 FOCUSS 并观察迭代过程

调用刚才写的 focuss 函数,参数取 p=0.5,max_iter=50,tol=1e-6:

x_hat = focuss(A, y, 0.5, 50, 1e-6); recovery_error = norm(x_hat - x_true) / norm(x_true); fprintf('相对重构误差: %.4f\n', recovery_error);

跑完以后,先把重构信号画出来,和真实信号叠加对比。正常情况下,你会看到 10 个真实非零位置处都有尖峰,而其他位置上的值接近 0,但不会严格等于 0——FOCUSS 不像硬阈值算法那样直接置零,它的稀疏性是"软"的。要硬约束的话,可以对结果做一步后处理:对 x_hat 取绝对值排序,保留前 K 个最大项,其余置零。这一步不在原始算法里,但很多工程场景会加。

我实测下来,同样的参数下 FOCUSS 的重构误差通常在 1e-4 量级。如果误差很大,优先检查是不是 p 值太小导致迭代发散,或者 A 没有归一化导致条件数过大。

3.3 把三种算法放到同一个平台上对比

资源包既然叫"几种常见算法",不对比就浪费了。我比较常放进来对比的是 OMP、BP 和 FOCUSS。OMP 在包里的实现一般长这样:

function x_hat = omp(A, y, K) [M, N] = size(A); r = y; support = []; x_hat = zeros(N, 1); for iter = 1:K % 计算残差与所有列的相关系数 corr = abs(A' * r); [~, idx] = max(corr); support = union(support, idx); % 最小二乘投影 x_temp = zeros(N, 1); x_temp(support) = pinv(A(:, support)) * y; r = y - A * x_temp; end x_hat = x_temp; end

BP 的实现稍微麻烦些,需要用线性规划或专门的 l1 求解器(比如 CVX 或 l1-MAGIC 工具箱)。如果资源包里没有相关工具箱,一个替代方案是用 CVX 的几行代码完成:

cvx_begin variable x_bp(N) minimize(norm(x_bp, 1)) subject to A * x_bp == y; cvx_end

对比时保持同一组 A 和 y,分别计算三种算法的重构误差和运行时间,结果整理成表格。典型结果会是:OMP 最快但依赖 K 先验,FOCUSS 居中且不需要 K,BP 误差最低但耗时最大。资源包里如果没带对比脚本,我建议你自己写一个,这是消化算法最好的方式。

3.4 加噪声场景下的参数调整

实际工程里观测一定有噪声。把观测改成 y = A*x_true + noise,其中噪声方差按信噪比 SNR 来设置。比如 SNR=20dB,对应噪声方差的公式是:

sigma_noise = norm(A * x_true) / sqrt(M) / (10^(SNR/20))

这在 MATLAB 里写着就是:

SNR_dB = 20; noise = randn(M, 1); noise = noise / norm(noise) * norm(A * x_true) / (10^(SNR_dB/20)); y_noisy = A * x_true + noise;

含噪时 FOCUSS 有个棘手问题:如果 p 太小,算法会把噪声也"稀疏化"一部分,导致解中出现虚假的非零元素。解决办法是把正则化微扰项调大,从 1e-6 提到 1e-3 甚至 1e-2。代价是重构精度下降,但解的稳定性明显增强。这本质上是稀疏性和抗噪性的权衡,没有绝对最优,只能根据实际信号的信噪比和稀疏度去调。

4. 常见问题与排查技巧实录

4.1 迭代发散结果全是 NaN 或 Inf

这个是我见过最多的问题。FOCUSS 的迭代式里有权重矩阵 W,而 W 的对角元是 |x|^(1-p/2)。如果某次迭代中 x 的某个分量恰好为 0,而 p<2,那么权重就会变成 0 的 1-p/2 次方——当 1-p/2 > 0 时结果是 0,没问题;但某些实现里如果出现负幂次,就会产生无穷大。解决方式有两个:一个是给权重加一个下界 epsilon,w = max(abs(x), eps).^(1-p/2);另一个是在更新 x 时加正则化项。资源包里的代码如果没做这个保护,建议自己加上。

排查时先用小规模数据测试,比如 N=64、M=16、K=3,逐行打印每次迭代的 x 范数,看是从哪一步开始异常的。通常问题都出在初始化阶段——如果最小二乘解 x_0 里有接近 0 的分量,第一轮迭代就容易踩坑。

4.2 重构误差很大但曲线形态正确

信号的大致轮廓画出来了,但具体数值对不齐。这种情况一般是两个原因:一是 p 值偏大导致稀疏性不够,很多本应置零的位置残留了小幅值;二是迭代次数不够,还没收敛就停了。先调大 max_iter 到 100 试试,如果误差明显下降,说明就是没跑够;如果误差没变化,再去调 p。

还有一种隐蔽情况:测量矩阵 A 的列没有归一化。如果 A 的各列范数差异很大,FOCUSS 会倾向于在大范数列方向集中能量,导致支撑集选错。检测方法很简单:计算 A 每列的范数,看是否接近相同值。资源包里的测试代码一般不会犯这种错,但如果你换成自己的测量矩阵,就一定要查。

4.3 OMP 效果好而 FOCUSS 差,怎么判断该用谁

这不是 bug,是特性。FOCUSS 本质上是在逼近 lp 范数解,它对测量矩阵的要求比 OMP 更严格。如果 A 的列之间相关性较强,OMP 靠逐步匹配反而更稳,FOCUSS 容易在多列之间摇摆。反之,如果 A 是严格随机的、列间相关性低的矩阵,FOCUSS 的重构精度通常优于 OMP,尤其是在稀疏度未知时。

实际项目里我的选择标准是这样的:如果稀疏度 K 已知或能估得很准,优先用 OMP,简单快;如果 K 未知且信号噪声适中,用 FOCUSS;如果追求极致精度且不心疼计算时间,用 BP/CVX。资源包把三种算法放一起,价值就在这儿——你可以在自己的数据上跑一遍对比,用数据说话,而不是凭感觉选。

4.4 速查表

问题现象可能原因排查方法解决建议
结果全是 NaN权重矩阵出现 0 的负次幂检查 W 对角元w = max(abs(x), eps).^(1-p/2)
重构误差大p 值过大或迭代不足调大 max_iter 观察减小 p,或增加迭代次数
误差差但形状对支撑集偏移检查 A 各列范数对 A 做列归一化
含噪时出现虚假尖峰p 过小导致噪声稀疏化检查解中非零分布增大正则化微扰项到 1e-3
运行极慢每次迭代显式求逆检查代码是否用 inv()改用反斜杠运算符或预分解
与 OMP 差距悬殊测量矩阵列相关性高计算列相关系数换用 OMP 或改门结构测量矩阵

4.5 调试技巧:在 FOCUSS 里埋观察点

调试这个算法有个小技巧:在每次迭代里临时存下 x 的非零位置和幅值,看它们随迭代的变化轨迹。如果算法正常,前几次迭代会有很多小幅值分量,随着迭代次数增加,这些分量逐步向接近 0 收缩,而少数主要分量保持增长或稳定。如果出现某个次要分量在后续迭代中反超主要分量,那就说明测量矩阵条件不好,或者 p 值出了问题。

我习惯的做法是每 5 次迭代输出一次当前稀疏度(即大于某个阈值的元素个数),观察它是否单调下降。FOCUSS 的稀疏度整体趋势是下降的,但中间会有波动,最终稳定在一个值附近。如果稀疏度不减反增,果断停掉调参,不要让它跑到最后。

5. 关于这个资源包,我的一点使用心得

从这个压缩包的命名就能看出来,它应该是一个实验性项目留下的产物——FOCUSS 是主菜,其他算法是配菜,全部打包在一起方便复用。我接触过不少这类资源包,最大的感受是:网上流传的代码质量参差不齐,但 FOCUSS 这类经典算法反而还算稳定,毕竟论文公开发表了三十年,核心公式不会有错,容易出问题的都在边界处理上。

我自己在实操中的体会是:拿到这类代码不要只跑 demo,一定要自己动手改参数、换数据、加噪声。FOCUSS 这个算法特别适合做这种实验,因为它就一个迭代式,参数就两三个,改动效果肉眼可见。你花一下午把 p 从 1 调到 0.2,跑一遍对比图,对稀疏重构的理解会超过干看一星期论文。最后再分享一个小技巧:调试稀疏重构算法时,先用小维度数据把所有参数摸清,再把维度拉到实际规模,否则定位问题会非常痛苦。希望这篇能帮你把这个压缩包彻底消化掉。

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

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询