企业官网建设流程全解析

简介：这是一份面向压缩感知（CS）学习与研究者的MATLAB代码包，聚焦多正弦信号的随机欠采样与重构问题。代码包内实现了正交匹配追踪（OMP）与SPGL1两种经典恢复算法，并配有示例脚本，能直观展示从稀疏建模、随机采样到迭代求解的完整链路，适合信号处理、通信工程等领域的初学者快速上手，也可供研究者在低采样率场景下对比算法性能。压缩包共36个文件，以25个m文件为主，涵盖算法主程序与示例函数；另含少量C/H源文件及预编译的mex文件，便于扩展和加速计算，整体仅67KB，轻量实用。目前已有1719人学习下载，代码目录结构清晰，运行示例即可理解压缩感知核心原理，为后续在数据采集、医学成像等方向的应用提供可直接修改的算法基底。

1. 压缩感知不是玄学：少采一半数据，照样把正弦信号找回来

拿到一列采样不完整的信号，缺失点位随机分布，直接做 FFT 会得到一坨乱谱——这是很多做频谱分析的人遇到过的头疼场景。压缩感知（Compressed Sensing, CS）正是在这种欠采样条件下，把背后真正的频率成分找出来的一套理论加算法组合。标题里写得很直白：多个正弦信号，随机欠采样，再用压缩感知恢复。核心价值很实际：采样点可以少于奈奎斯特率要求，只要信号在某个基下是稀疏的——正弦信号天然在傅里叶基下稀疏，所以这个问题是压缩感知最典型的应用入口。做信号处理实验的、搭数据采集链路的、做振动监测或频谱监测的工程师，都能在这套流程里找到可直接复用的 Matlab 实现。

本文不堆公式，直接从“测量矩阵怎么构造”讲到“OMP 怎么实现”，再讲到稀疏度、字典网格、幅度相位归一化这些一调就翻车的细节。代码按小节给全，复制到 Matlab 里能跑通，参数也给出经验取值和调整方向。

2. 随机欠采样怎么设计：测量矩阵与采样 pattern 的取舍

压缩感知的恢复质量，有一半在采样端就决定了。正弦信号的恢复，测量对象是频谱系数，但采样动作发生在时域——这种“时域随机抽点、频域看稀疏”的结构，决定了测量矩阵必须选对。

2.1 为什么用傅里叶测量矩阵而不是高斯随机矩阵

理论上随机高斯矩阵的约束等距性（RIP）最好，对稀疏信号几乎“什么都能恢复”。但在正弦信号这个具体问题上，我更倾向于用部分傅里叶矩阵，也就是“从完整 DFT 矩阵里随机抽取若干行”构成的测量矩阵。原因有两个：一是物理意义直接，时域随机欠采样点经过 DFT 行抽取，等价于观测了频谱的部分信息，硬件上就是把采样时刻打乱，不需要额外做矩阵乘法；二是正弦信号在傅里叶基下本来就稀疏，测量基和稀疏基天然配对，恢复时相关运算用 FFT 可以加速，不必显式构造大矩阵。

下表是对三种常见测量矩阵的选型判断。

我实际做正弦恢复时，默认用部分傅里叶矩阵。高斯矩阵只有在换稀疏基做验证时才会用，比如时域脉冲信号的恢复。

2.2 随机索引的两种生成方式：完全随机与分块随机

采样 pattern 直接影响恢复成功率。最简单的做法是randperm随机抽 m 个时域样本，这种方式对理论推导友好，但实际中可能出现两个问题：某一段连续样本全部缺失，形成长空窗；或者某一段样本过密，信息冗余。工程上我更推荐分块随机：把 N 个样本分成若干块，每块内部随机抽固定数量，既保证随机性，又避免长空窗。

N = 512; % 采样总点数 m = 128; % 欠采样点数，约为总点数的 25% rng(2024); % 固定随机种子，便于复现 % 方式一：完全随机索引 meas_idx_full = sort(randperm(N, m)); % 方式二：分块随机索引 nb = 16; % 分块数，经验值取 N/32 左右 blockLen = N / nb; % 每块长度 r = ceil(m / nb); % 每块抽样数，m=128 时每块 8 个 meas_idx_block = zeros(1, r * nb); for b = 1:nb offset = (b - 1) * blockLen; pick = offset + randperm(blockLen, r); meas_idx_block((b - 1) * r + 1 : b * r) = pick; end meas_idx_block = sort(meas_idx_block(1:m)); % 裁剪到 m 个点

这段代码的关键在分块参数：nb选 16 时，每块 32 个点，块内抽 8 个，块与块之间不会出现整块全部缺失的情况。randperm(blockLen, r)保证块内索引不重复，sort只是让索引有序，方便画图观察 pattern，对后续构造矩阵没有影响。完全随机方式适合理论对比，分块随机方式适合实际采样链路，尤其是那些一次采集时间窗很短、不能容忍长空窗的场景。

2.3 从欠采样到压缩感知问题的数学形式

把信号 x 视为 N 维时域向量，x 在傅里叶基下可以表示为 x = Ψs，其中 Ψ 是归一化 IDFT 矩阵，s 是频域稀疏系数。随机欠采样的操作是保留部分时域位置的值，等价于在等式两边同时乘以一个行抽取矩阵 R，得到观测模型：

y = RΨs = A s

这里 y 是 m 维观测向量，A = RΨ 是 m×N 的测量矩阵，也就是部分傅里叶矩阵。压缩感知恢复要做的就是：已知 A 和 y，在 s 稀疏的约束下解出 s。这个模型里有两个容易忽略的细节：一是 Ψ 必须归一化，否则幅度恢复会差一个比例因子；二是 R 固定后 A 就确定了，OMP 每次迭代都要在 A 的列上做内积，所以 A 的存储和计算效率直接影响恢复速度。

3. 用 Matlab 把恢复跑通：OMP 实现与幅度相位提取

压缩感知的恢复算法很多，基追踪、ISTA、FISTA、OMP 各有适用面。正弦信号恢复这种稀疏度明确的问题，我首选 OMP，因为它直观、迭代次数可控、参数少，对频谱泄漏不敏感的场景非常稳定。

3.1 OMP 主循环：相关、选原子、最小二乘、更新残差

OMP 的核心逻辑是贪心：每次迭代找与当前残差相关性最强的原子，把它加入支撑集，然后在支撑集上重新最小二乘，更新残差，直到迭代次数用尽或残差足够小。

function s_hat = omp(A, y, K, tol) % A: m*N 复数测量矩阵 % y: m*1 观测向量 % K: 稀疏度（预期非零谱线数） % tol: 残差能量停止阈值，传 0 表示不启用，跑满 K 次 [m, N] = size(A); s_hat = zeros(N, 1); r = y; % 残差初始化为观测向量 support = []; % 支撑集，存放已选原子索引 x_temp = []; for iter = 1:K % 第 1 步：残差与所有原子的相关 corr = A' * r; % N 维相关向量，A' 是共轭转置 [~, idx] = max(abs(corr)); % 选相关最强的原子 support = union(support, idx); % 加入支撑集，union 防重复 % 第 2 步：在支撑集上对原观测 y 做最小二乘 At = A(:, support); x_temp = At \ y; % 第 3 步：用新系数更新残差 r = y - At * x_temp; % 第 4 步：提前判断收敛 if tol > 0 && norm(r) < tol break; end end s_hat(support) = x_temp; end

这里最关键的一步是第 2 步的At \ y——用原观测 y 而不是残差 r 做最小二乘。很多人第一次写 OMP 会用At \ r，结果每次迭代的系数是局部修正量，支撑集一旦包含两个相关原子就乱套。OMP 的正交性来自“每次选完原子后重新对 y 投影”，不是对残差投影。tol参数在有噪声场景下建议设成噪声能量的量级，比如1e-3；纯信号场景设1e-6以上，避免低信噪比时迭代过头。

3.2 主程序：三个正弦信号从欠采样到恢复

下面这段主程序把采样、欠采样、恢复、画图串起来，全部代码量不超过 60 行，在普通笔记本上跑完不到一秒。

N = 512; % DFT 点数 Fs = 1024; % 采样率，Hz t = (0:N-1) / Fs; % 三个正弦：幅度分别 1.0 / 0.6 / 0.3，频率 37 / 157 / 211 Hz x = 1.0 * sin(2*pi*37*t) + 0.6 * sin(2*pi*157*t) + 0.3 * sin(2*pi*211*t); rng(7); x = x + 1e-3 * randn(size(x)); % 加一点噪声，避免零残差假象 K = 3; % 稀疏度：三个正弦 m = 128; % 欠采样样本数，采样率只有 12.5% rng(2024); meas_idx = sort(randperm(N, m)); % 随机欠采样索引 % 构造部分傅里叶测量矩阵，注意除 sqrt(N) 归一化 Psi = exp(1j * 2 * pi * meas_idx(:) * (0:N-1) / N) / sqrt(N); y = x(meas_idx); % 时域观测：只保留采样点 % 恢复频谱系数 s_hat = omp(Psi, y, K, 1e-6); % 画图对比 freqs = (0:N-1) / N * Fs; stem(freqs, abs(s_hat), 'filled'); xlabel('Hz'); ylabel('幅度');

这段代码里的测量矩阵构造是重点：meas_idx(:)是 m×1 列向量，(0:N-1)是 1×N 行向量，两者相乘得到一个 m×N 的指数矩阵，除以sqrt(N)后保证能量归一化。如果不除sqrt(N)，恢复出的幅度整体偏大 N 倍，后面幅度提取全部作废。K=3告诉 OMP 只找三根谱线，这个前提在已知信号结构时合理；如果真实信号里还有直流或未知谐波，K 要相应加大，第 4 章会讲怎么估计。

3.3 从恢复系数反推幅度与相位：归一化别搞错

OMP 恢复出来的是复频谱系数，不是直接的时域幅度。取模再除以归一化因子才是正弦幅度，这一点几乎每个第一次做的人都会踩。用上面sqrt(N)归一化的测量矩阵时，正交 DFT 使正弦幅度与频谱系数满足：

幅度 = 2 × |s_peak| / sqrt(N)

直流分量不乘 2，直接是 |s_peak| / sqrt(N)。相位方面，正弦信号 sin(2πft + φ) 的频谱系数，其辐角等于 φ - π/2，所以相位要加回 π/2。

% 提取峰值索引，按幅度降序 [amp_sorted, peak_idx] = sort(abs(s_hat), 'descend'); for p = 1:length(peak_idx) idx = peak_idx(p); if abs(s_hat(idx)) < 1e-3 * max(abs(s_hat)) % 忽略数值噪声 continue; end if idx == 1 % 直流分量 amp = abs(s_hat(idx)) / sqrt(N); else amp = 2 * abs(s_hat(idx)) / sqrt(N); end phase = angle(s_hat(idx)) + pi/2; % 对 sin 形式的相位修正 if phase > pi, phase = phase - 2*pi; end fprintf('f=%.1f Hz, amp=%.3f, phase=%.3f\n', freqs(idx), amp, phase); end

这个提取逻辑依赖于 DFT 定义的一致性：测量矩阵用sqrt(N)归一化，频谱系数就是正交 DFT，幅度公式才成立。如果你习惯用 Matlab 的fft不归一化，那幅度公式要改成2 * abs(s_hat(idx)) / N，相位修正不变。我建议固定用sqrt(N)归一化这一套，写进自己的工具函数里，省得每次换算出错。

4. 稀疏度、字典网格与正则化：三个必调参数

OMP 表面参数少，但稀疏度 K、字典网格密度、是否引入正则化这三件事，决定了恢复结果是从“能看”到“能用”的分水岭。

4.1 稀疏度 K 怎么估：残差膝盖点法

K 是 OMP 停止迭代的直接条件。K 设小了，弱分量直接被漏掉；K 设大了，后面几次迭代在纯噪声上找原子，恢复结果出现虚假峰。已知信号结构时按真实谱线数设 K 即可；不知道结构时，我常用“残差膝盖点”来定 K：跑若干次迭代，记录每次残差能量，画semilogy曲线，曲线从陡降转为平缓的拐点就是 K。

maxIter = 30; res = zeros(maxIter, 1); for k = 1:maxIter s_tmp = omp(Psi, y, k, 0); % tol 传 0，强制迭代 k 次 res(k) = norm(y - Psi * s_tmp); end figure; semilogy(1:maxIter, res, 'o-'); xlabel('迭代次数'); ylabel('残差能量'); grid on;

运行这段脚本会看到：K=3 之前残差指数级下降，K=3 之后残差降到噪声底并基本不再变化，膝盖点就是真实稀疏度。配合实际场景，如果信号里有微弱谐波淹没在噪声里，膝盖点会变得模糊，这时要结合频谱背景噪声水平决定 K，不要追求残差无限小，否则恢复的是噪声。

4.2 字典网格密度：频率分辨率与 RIP 的拉锯

OMP 默认在 N 点 DFT 网格上找频率，频率分辨率是 Fs/N。如果信号频率正好落在网格上，恢复漂亮；如果不落在网格上，稀疏假设被破坏，一根正弦会泄漏到多个网格原子。常见做法是加密字典网格，比如把网格数从 N 扩到 N2 = γN，频率分辨率提高到 Fs/(γN)。

gamma = 4; % 网格加密倍数 N2 = gamma * N; freq_grid = (0:N2-1) / N2 * Fs; Psi2 = exp(1j * 2 * pi * meas_idx(:) * (0:N2-1) / N2) / sqrt(N); % 恢复时 OMP 的支撑索引对应 freq_grid s_hat2 = omp(Psi2, y, K, 1e-6);

网格加密有代价：字典原子之间的相关性变高，测量矩阵 RIP 性质变差，恢复所需的采样点数 m 必须相应增加。经验上 m 与稀疏度 K、字典规模 N2 满足 m ≈ 2~4 × K × log(N2/K)，γ 超过 8 后提升有限，计算量却线性上涨。实际工程里我会先用 γ=1 快速跑一遍，观察峰值是否落在网格边界上，必要时只对峰值邻域做局部加密，而不是全字典加密。

4.3 OMP 为什么不需要正则化参数，基追踪才需要

OMP 的每一步都在支撑集上做正交投影，系数由最小二乘唯一确定，所以没有正则化参数。但基追踪（BP / LASSO）在稀疏约束和拟合误差之间加了一个 λ，λ 的值对结果影响很明显，需要做交叉验证或按理论公式估计：λ ≈ 噪声标准差 × sqrt(2 × log(N))。在正弦恢复场景，如果噪声水平很低、稀疏度明确，OMP 是更省心的选择；如果频带内有大量弱分量、稀疏度不好定，基追踪能一次给出分布式的解，不会把能量集中在少数峰值上。我的习惯是先用 OMP 摸清结构，再用基追踪验证结果，两者互相印证比单用任意一个都靠谱。

5. 避坑指南：正弦恢复最容易翻车的四个现场

压缩感知恢复流程不长，但每一步都有一个“看着没问题、跑出来全错”的坑。以下四条是我实际调试中反复踩过的，每条按现象、原因、解决写清楚。

5.1 低频与直流偏置把稀疏度名额占掉

现象：信号里有个直流偏置，K 设成 3，恢复结果把频率最高的弱分量丢了，换回来一根不在任何正弦附近的谱线。原因：直流在频谱里能量占比很高，OMP 第一次迭代就选它，3 次迭代里两次被 DC 和主强分量占掉，弱分量根本没机会。解决：恢复前先把观测信号的直流均值减掉，把x = x - mean(x)放在采样前；如果底层传感器确实有偏置，把 K 加 1，并在支撑集输出里手动忽略 DC 索引。

5.2 两个频率挨太近，OMP 把原子选错

现象：两个正弦频率只差一个 FFT bin 多一点，恢复出的两根谱线和真实频率对不上。原因：字典里相邻原子的相关性太高，残差对两个原子的投影区分不明显，贪心算法在第一步就可能选错。解决：缩短采样时长无法改善频率分辨率时，先做粗网格恢复确定大致频带，再在粗峰值邻域用更细网格做局部搜索，相当于两段式估计；同时可以尝试把频率差小于 2×Fs/N 的两根正弦不再区分，而是合并为一个宽带模型再拟合。

5.3 完全随机抽样出现长空窗，低频线索丢失

现象：每次运行rng不同种子，恢复结果时好时坏，低频分量偶尔大幅偏差。原因：完全随机索引在样本量小时，可能出现一段连续几十个点全部缺失。低频正弦变化缓慢，本来就靠长时间累积体现，空窗期正好把这段信息丢光。解决：采用第 2.2 节的分块随机 pattern，把长空窗从机制上排除；采样点数 m 紧贴理论下限时，多跑几次不同种子，取恢复成功率而不是单次结果下结论。

5.4 频率不在网格上，恢复系数“一峰变三峰”

现象：信号是 157.3 Hz，网格分辨率 2 Hz，恢复结果在 156 Hz、158 Hz 附近出现两个小峰，主峰幅度只有真实值的三分之二。原因：正交字典假设频率恰好在网格点上，off-grid 频率无法被单个网格原子稀疏表示，能量泄漏到相邻原子。解决：先用网格恢复确定主峰位置，再在主峰两侧做抛物线插值估计真实峰位；更彻底的做法是改用连续字典或网格自适应迭代，但这会让恢复时间上升一个量级，工程上优先用插值修正。幅度修正公式是真实幅度 ≈ 测量峰值幅度 / sinc(δ)，δ 是真实峰位与网格点的归一化偏差。

6. 验证不是看一眼波形图：用蒙特卡洛曲线判断能不能上线

单次恢复成功说服不了任何人，尤其当采样点数紧贴理论极限时，压缩感知的成功率本身就是个随机变量。我判断这套恢复流程能不能用，不看恢复出的波形漂不漂亮，而是看蒙特卡洛成功率曲线。

6.1 从单次恢复误差到批量成功率

mSet = 64:32:320; % 不同欠采样点数 numTrial = 50; % 每个 m 下的随机试验次数 succRate = zeros(size(mSet)); for mi = 1:length(mSet) m = mSet(mi); succCnt = 0; for trial = 1:numTrial % 随机生成三根正弦，频率在网格上随机抽取但互不邻近 f_true = sort(randperm(200, K) + 50); x = zeros(1, N); for k = 1:K x = x + (0.5 + 0.4*rand) * sin(2*pi*f_true(k)*t); end meas_idx = sort(randperm(N, m)); Psi = exp(1j * 2 * pi * meas_idx(:) * (0:N-1) / N) / sqrt(N); s_hat = omp(Psi, x(meas_idx), K, 1e-6); [~, idx] = sort(abs(s_hat), 'descend'); est = sort(idx(1:K)); if isequal(est, f_true) % 频率索引完全一致才判成功 succCnt = succCnt + 1; end end succRate(mi) = succCnt / numTrial; end figure; plot(mSet, succRate, 'o-'); ylabel('恢复成功率'); xlabel('m');

判成功标准用“频率索引完全一致”是很苛刻的，但只要网格点足够细、噪声足够低，它最能暴露参数问题。跑出来的曲线通常是三段：m 很小时成功率在 20% 以下，中间一段快速爬升，超过某个阈值后稳定在 100%。这个阈值 m_min 就是你的系统在这个稀疏度下能可靠工作的最小采样点数。

6.2 成功率曲线怎么解读

曲线拐点位置比单点成功率有用。假设 N=512、K=3，理论下界约是 m ≈ 3K·log(N/K) ≈ 47，实际曲线拐点通常在 100 到 160 之间，差距来自随机测量矩阵的方差和网格相关性。如果你的应用场景能接受 95% 成功率，就选拐点右侧的 m；如果系统有实时采样率上限，则反过来看当前 m 能达到多高的成功率。

我自己的习惯是把这条曲线画成默认实验的一部分，换信号模型、换字典、换测量 pattern 时都会重跑一遍。曲线的形状比那条恢复出的正弦波形更诚实，因为它把随机种子、低信噪比这些不可控因素全部摊开来看，而不是靠一次好运气。这套流程从采样 pattern 到 OMP 再到验证曲线，全部代码不超过两百行，改几个参数就能用于你自己的采集链路。希望帮到你。

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