MATLAB压缩感知实战：随机欠采样与OMP恢复

简介压缩感知(Compressed Sensing, CS)的Matlab实现代码包专注多正弦信号的随机欠采样与精确重构面向信号处理研究者与工程技术人员适合想快速验证CS理论的学习者。包内基于稀疏表示框架提供正交匹配追踪(OMP)与SPGL1两种恢复算法通过随机采样策略突破奈奎斯特限制并配有演示脚本便于对比两种算法在噪声环境与低过采样率下的重构表现。资源共36个文件以25个m源程序为主体辅以C/H源码和MEX编译文件整体67KB体积精简、目录清晰可直接在Matlab中运行。已有1709人学习/下载可用于采样矩阵设计、重构算法评估也为数据采集、无线通信、医学成像等场景提供了实践参考。1. 从随机欠采样说起压缩感知为什么能救回丢失的正弦信号假设你手里只有一个多正弦叠加信号在 512 个均匀时刻的采样值但采集过程中受存储或硬件限制实际只保留了其中随机抽出的 128 个时刻的幅度其余 384 个点全部丢失。直觉上这是一个“缺了一大半数据”的插值问题但直接插值会完全失败。压缩感知给出了一条反直觉路径只要信号在某个变换域比如傅里叶频域足够稀疏随机欠采样得到的数据其实包含恢复所需的全部信息真正要做的不是补点而是求解一个稀疏系数向量。这个结论直接改变了雷达、超声成像和高频宽带采样系统的设计思路。接下来我用 MATLAB 把多正弦信号的随机欠采样、观测矩阵构造、OMP 恢复整体走一遍并标出参数设置的边界和常见误区。2. 压缩感知的建模前提稀疏字典与随机观测矩阵压缩感知的第一步是把“时域丢点”这个现象改写成一个线性观测模型。设原始等间隔采样的信号为长度为 N 的列向量 $x$观测矩阵 $\Phi$ 是一个 $M \times N$ 的行选择矩阵$M \ll N$实际观测向量 $y \Phi x$。如果 $x$ 本身不稀疏就需要找到一个正交变换矩阵 $\Psi$使得 $x \Psi s$其中 $s$ 只有 $K$ 个非零元素那么观测方程就变成 $y \Phi \Psi s$。2.1 多正弦信号的傅里叶字典用 DFT 矩阵还是显式构造多正弦信号在频域天然是稀疏的因此 $\Psi$ 首选 DFT 矩阵。MATLAB 中有两种常见做法。一种做法是直接调用dftmtx(N)得到完整的 $N \times N$ 傅里叶变换矩阵再取被采样行。dftmtx返回的矩阵 $F$ 满足 $FF N I$逆变换要用 $F$ 除以 $N$。另一种做法是显式构造正弦/余弦原子库把字典列设为不同频率下的 $\cos(2\pi f t)$ 和 $\sin(2\pi f t)$。两种方式的差别在于字典类型稀疏系数字典大小典型场景DFT 复指数字典复数且正负频率成对出现$N \times N$部分傅里叶观测CS 理论中最常用的模型实数正弦/余弦字典实数每个频率占两列$2F \times N$需要强制实系数、且不希望出现共轭对时我通常用 DFT 矩阵做验证实验因为它与 FFT 完全对应便于用fft(x)的结果对比恢复出的系数。设 $F$ 为 DFT 矩阵则 $s Fx$ 是频域系数。对于三个正弦信号$s$ 中只有 6 个非零主峰每个频率对应正负两条谱线符合稀疏条件。2.2 随机欠采样矩阵从随机下标到部分傅里叶观测矩阵随机欠采样的操作非常简单在 $1$ 到 $N$ 中随机抽取 $M$ 个下标作为保留位置。对应的选择矩阵 $\Phi$ 是一个稀疏矩阵每一行只有一个 1其余为 0第 $i$ 行记录第 $idx(i)$ 个采样点在原始信号中的位置。把选择矩阵与 DFT 矩阵相乘得到一个部分傅里叶矩阵 $A \Phi F$。此时观测向量为$$ y A s $$其中 $s$ 是稀疏频域向量。实际代码里不需要显式生成稀疏的 $\Phi$直接取 $F$ 的若干行即可。随机选择下标的意义在于将频域采样引起的混叠转化为类似噪声的干扰。如果按等间隔丢掉数据丢失的频点会产生严重的频谱泄漏难以用稀疏优化恢复随机丢弃时重叠的旁瓣被随机化压缩感知的恢复算法可以把真正的稀疏分量从“噪声背景”中挑出来。2.3 为什么用 OMP 而不是直接解最小二乘观测方程 $y A s$ 是一个欠定方程组未知量 $N$ 个方程只有 $M$ 个直接最小二乘会得到能量分散的非稀疏解。压缩感知的核心是用稀疏性先验来约束解典型思路是最小化 $\ell_0$ 范数但这是个 NP-hard 组合问题。实际实现中使用 $\ell_1$ 范数凸松弛或者用迭代贪婪算法近似求解。优化形式优化目标计算复杂度适用条件$\ell_0$$\min |s|_0$ s.t. $Asy$NP-hard只适合理论分析$\ell_1$$\min |s|_1$ s.t. $Asy$凸优化可解需要 CVX、yall1 等工具OMP迭代选择相关性最大的原子约 $O(KMN)$稀疏度已知或可估计工程中最常用正交匹配追踪OMP的思想是每次迭代从 $A$ 的所有列中找与当前残差相关性最强的一列加入支撑集然后用最小二乘更新支撑集上的系数重新计算残差重复 $K$ 次。这样避免了直接求 $\ell_1$ 凸优化时对第三方优化工具箱的依赖也更容易在单片机或 DSP 上实现。3. 在 MATLAB 中实现多正弦信号的随机欠采样与 OMP 恢复现在把建模过程变成可运行的 MATLAB 代码。下面的示例在N 512个均匀采样点上生成三个正弦信号的叠加随机抽取 $M 128$ 个点然后通过 OMP 从部分傅里叶观测中恢复频域系数和原始时域波形。3.1 生成多正弦信号并设置频率格点定义采样率fs和总点数N让每个正弦频率都落在 DFT 频率分辨率的整数倍上否则会引入频谱泄漏增大恢复难度。N 512; % 原始等间隔采样点数 fs 1024; % 奈奎斯特采样率 t (0:N-1) / fs; % 时间列向量N x 1 f0 [50, 120, 200]; % 三个正弦频率单位 Hz A0 [1.0, 0.8, 0.6]; % 对应幅度 % 频率分辨率 fs / N 2 Hz所有频率都被 2 整除位于FFT网格上 x A0(1)*sin(2*pi*f0(1)*t) ... A0(2)*sin(2*pi*f0(2)*t) ... A0(3)*sin(2*pi*f0(3)*t);代码中把x设计为列向量与后续dftmtx的矩阵乘法维度保持一致。三个频率对应频域中的 6 条谱线所以真正的稀疏度是 6。如果信号的频率不是频率分辨率的整数倍比如f0 [50.5, 120, 200]那么频谱不再稀疏OMP 的恢复误差会显著增大这一点后面会专门说明。3.2 随机欠采样与部分傅里叶观测矩阵随机欠采样的关键是使用randperm生成不重复的下标。为了实验可复现我习惯用rng固定随机种子。M 128; % 实际保留的采样点数 rng(42); % 固定随机种子保证结果可复现 idx sort(randperm(N, M)); % 随机抽取 M 个不同采样下标并排序 y x(idx); % 实际观测值 F dftmtx(N); % 完整 N x N DFT 矩阵 A F(idx, :); % 部分傅里叶观测矩阵M x N这里的A就是前面说的 $A \Phi F$。dftmtx(N)生成的是复数矩阵因此A也是复数后续 OMP 中相关运算必须使用模值。观测值y是实数但被复矩阵投影后解空间需要按复数处理。OMP 找到的稀疏系数 $s$ 会是一个复向量其中正负频率各有一对共轭对称的峰值。3.3 OMP 函数与信号重建下面是一个可以直接放到脚本或独立.m文件中的 OMP 实现。function s_hat omp(A, y, K) % 正交匹配追踪从部分傅里叶观测 y 中恢复稀疏系数 s % A: M x N 观测矩阵 % y: M x 1 观测向量 % K: 稀疏度预期非零系数个数 M size(A, 1); N size(A, 2); r y; % 残差初始为观测向量 support []; % 支撑集记录已选原子的索引 s_hat zeros(N, 1); for iter 1:K % 排除已选原子计算所有剩余原子与残差的相关系数 remaining setdiff(1:N, support); corr A(:, remaining) * r; % 每个原子与残差的内积 [~, j] max(abs(corr)); % 取模值最大者j 在 remaining 中的位置 idx_j remaining(j); % 转换为全局原子索引 support [support, idx_j]; % 加入支撑集 % 用最小二乘重新估计支撑集上所有原子的系数 A_s A(:, support); coef A_s \ y; % 求解超定方程的最小二乘解 r y - A_s * coef; % 更新残差 if norm(r) 1e-10 % 残差足够小提前停止 break; end end s_hat(support) coef; endremaining setdiff(1:N, support)保证了同一原子不会被重复选入避免支撑集大小在迭代中无意义增长。corr的长度等于剩余原子数量而remaining(j)将局部索引映射回全局字典列索引。A_s \ y在 MATLAB 中会自己选择合适的最小二乘算法因为A_s是典型的小规模矩阵速度足够快。恢复频域系数后用逆 DFT 重建时域信号并计算相对误差K 6; % 三个正弦对应 6 条谱线稀疏度 6 s_hat omp(A, y, K); x_hat ifft(s_hat, N); % 因为 s_hat 是 F*x 的近似ifft 即逆变换 err norm(x - x_hat) / norm(x); % 相对误差 fprintf(相对重建误差: %.4e\n, err);ifft默认对向量作归一化逆变换与dftmtx的定义一致。注意这里使用的是ifft而不是A或F因为dftmtx生成的是未归一化 DFT 矩阵而ifft自带 $1/N$ 缩放因子正好得到时域信号。如果恢复成功x_hat和原始x在时域几乎重合相对误差通常在1e-6以下如果M过小或者K估计错误误差会骤然升高。4. 参数边界稀疏度 K、采样点数 M 和字典规模怎么搭配CS 并不是无条件成立的。即使信号本身稀疏观测矩阵也需要满足一定约束等距性质RIP才能保证恢复。对于部分随机傅里叶矩阵取得的理论条件是 $M \geq C K \log(N/K)$其中常数 $C$ 一般在 2 ~ 4 之间。实际工程中需要根据实验扫描确定参数下限。4.1 用一组经验参数规避理论死角RIP 条件很严格实际调试时更常用的是“采样点数 vs 稀疏度”的经验比例。对于随机部分傅里叶观测矩阵以下参数范围可作为起点信号长度 N稀疏度 K随机采样点数 M恢复效果预期512232稳定恢复误差 1e-6512448稳定恢复误差 1e-6512696多数随机种子可恢复边界区域5126128稳定恢复余量充足102410160稳定恢复但计算时间明显上升这张表不是理论下界而是我做过多次随机实验后得到的“安全区”。可以看到稀疏度从 2 升到 6要求采样点数远不是简单的 3 倍关系因为 $\log(N/K)$ 带来的非线性影响在低 $M$ 时非常明显。4.2 扫描 M/N 以定位恢复临界点用一段循环代码对不同 $M$ 进行测试是定位临界点最直接的方式。下面的脚本固定随机种子逐个扫描M观察相对误差的拐点。M_list [32, 48, 64, 80, 96, 112, 128]; % 不同保留点数 err_list zeros(size(M_list)); for i 1:length(M_list) M_i M_list(i); rng(42); idx_i sort(randperm(N, M_i)); A_i F(idx_i, :); y_i x(idx_i); s_i omp(A_i, y_i, K); x_i ifft(s_i, N); err_list(i) norm(x - x_i) / norm(x); end运行后观察err_list通常会出现一个明显的分界在 $M$ 较小时误差接近 1表示恢复完全失败增大到某个阈值后误差突然降到1e-6以下这个拐点就是当前信号和字典条件下的临界采样率。这里必须特别注意rng(42)每次只固定了随机采样下标OMP 本身的随机性取决于字典和信号因此结果可以复现。4.3 频率离格时如何处理如果正弦频率不是 FFT 格点的整数倍比如 50.5 Hz 而不是 50 HzFT 字典中没有任何一个原子能精确表示它。此时 $s$ 不再稀疏OMP 会把能量泄漏到相邻的若干频率单元上恢复误差很大。常见做法是缩小频率间隔把 $N$ 增大到 1024 或 2048相当于加密频域网格。但代价是 $A$ 变为 $M \times 2048$OMP 中每次内积运算的计算量随之增长。另一种做法是先做一次 FFT 粗估计频率位置然后在估计频率附近局部细化网格把全局稀疏问题转成若干个局部稀疏问题这种多分辨率思路在工程中很常用。5. 延伸与验证把压缩感知恢复做成可诊断的测量流程OMP 代码能跑通只是第一步实际使用中更关键的是确认恢复结果是否可信。一个简单而有效的校验方法是观察残差的能量与支撑集的变化。如果残差在迭代过程中下降到噪声水平后不再下降说明支撑集大小基本饱和如果残差在若干次迭代后仍然与初始观测向量同量级大概率是M过小或稀疏度K估计偏大。另一个诊断技巧是重复运行多次随机欠采样比较每次恢复出的频率位置是否一致如果频率支撑对随机种子极其敏感说明当前条件接近恢复失败边界需要增加采样点数。更进一步可以把静态压缩感知扩展成自适应压缩感知。静态方案一次性选好 $M$ 个采样点而自适应方案先取一小部分采样点用恢复结果估计信号能量集中在哪些频段再针对高能量区域增加局部采样。这样做的好处是在总采样率受限时能把宝贵的采样资源集中在真正有信号的子带。实现上并不复杂只需要在第一轮恢复后得到支撑集索引再优先采集靠近这些频率的时域点第二轮重新做 OMP 时可以把新老观测合并进同一个 $y$ 和同一个部分傅里叶矩阵中。对于慢变化的多正弦信号这种两阶段测量比均匀随机抽取的恢复成功率更高。实际调试时建议第一轮固定rng(42)把实验条件完全复现出来确认算法正确后再逐个松开随机种子观察统计意义上的恢复成功率。如果使用dftmtx生成的大矩阵导致内存紧张可以用A exp(2i*pi*(0:M-1)*(0:N-1)/N)直接构造部分傅里叶矩阵但必须注意频率索引与fft输出位置一致避免正频率和负频率的顺序错位。最后做 CS 恢复之前先画一次原始信号和欠采样点的散布图直观确认采样下标确实覆盖了完整时间范围如果随机种子设置不当导致采样点集中在前半段矩阵 $A$ 的各行强相关恢复就会失败。这一步检查经常能省下大量排错时间。本文还有配套的精品资源点击获取