企业官网建设流程全解析

简介：压缩感知(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$ 满足 $F'F = N I$，逆变换要用 $F'$ 除以 $N$。另一种做法是显式构造正弦/余弦原子库，把字典列设为不同频率下的 $\cos(2\pi f t)$ 和 $\sin(2\pi f t)$。两种方式的差别在于：

我通常用 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$ 范数凸松弛，或者用迭代贪婪算法近似求解。

正交匹配追踪（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; end

remaining = 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 稀疏度”的经验比例。对于随机部分傅里叶观测矩阵，以下参数范围可作为起点：

这张表不是理论下界，而是我做过多次随机实验后得到的“安全区”。可以看到，稀疏度从 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 Hz，FT 字典中没有任何一个原子能精确表示它。此时 $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$ 的各行强相关，恢复就会失败。这一步检查经常能省下大量排错时间。

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