企业官网建设流程全解析

简介：一份以压缩感知（Compressed Sensing）为核心的多正弦信号恢复MATLAB代码包，面向信号处理、通信、医学成像等领域的研究者与工程师，帮助理解并实践远低于奈奎斯特采样率的随机欠采样与稀疏重构方法。资源共36个文件、压缩包约67KB，以25个m脚本为主体，配合c与h源码、mexw32等动态链接文件，覆盖从基础演示到核心算法（正交匹配追踪OMP、SPGL1）的完整调用链路。已有1709人下载学习。包内代码不仅实现了多个正弦信号的随机欠采样与重构，还给出了L1范数最小化、L1/L2范数投影等辅助函数，便于对比OMP与SPGL1在不同噪声环境和稀疏度下的恢复效果。通过运行示例，读者可以直观掌握压缩感知的建模流程、参数调节思路以及重构质量评估方法，为一维稀疏信号的工程化采集和恢复提供可直接改造的参考脚本。

1. 压缩感知打破奈奎斯特采样定理的限制

ADC 采样率告急时，工程师的第一反应是换更高速的芯片，但很多时候算法可以替代硬件——压缩感知（Compressed Sensing, CS）恰恰证明了这一点。它让我们在远低于奈奎斯特率的条件下随机获取少量采样点，仍然能精确恢复原始信号，前提是信号在某个变换域足够稀疏。这个反直觉的结论不是建立在稀疏字典的花哨包装上，而是落在最普通的一个事实里：一个由几个正弦波叠加的信号，在频域里只是寥寥数个非零系数。你真正需要测量的，不是信号本身的全部时间样本，而是这些极少数非零频点在测量矩阵上的投影。本文将用 MATLAB 从零实现多正弦信号的随机欠采样与恢复，涉及随机测量矩阵构造、稀疏基选择、OMP 与基追踪算法的具体实现和参数调节，手把手跑通一整套可复现流程。

2. 压缩感知恢复正弦信号的三块基石：稀疏表示、随机测量与恢复算法

2.1 正弦信号在傅里叶基下的稀疏性为什么是恢复的前提

压缩感知的第一条铁律是信号必须可稀疏表示。对连续时间信号采样得到离散序列 x[n]，如果它由 K 个不同频率的正弦波叠加而成，那么对其做离散傅里叶变换（DFT），频谱只在 K 个频点处有非零值。

N = 512; % 信号长度 fs = 1000; % 采样率 1000Hz t = (0:N-1)/fs; % 时间轴 f = [50, 150, 267]; % 三个正弦频率 x = sin(2*pi*f(1)*t) + 0.8*sin(2*pi*f(2)*t) + 0.5*sin(2*pi*f(3)*t); X = fft(x)/N; % 归一化 FFT stem(abs(X(1:N/2))); % 只看单边谱

这段代码把三个频率的正弦波叠加成时域信号，再通过 FFT 转到频域。跑完后你会看到谱线只在 50、150、267Hz 三个位置凸起，其余位置接近机器精度。K=3 就是信号的稀疏度，也是后面测量矩阵设计、恢复算法迭代次数设定的核心依据。稀疏度估计偏低时恢复误差迅速恶化，偏高时计算量浪费但结果仍然正确，这与匹配追踪类算法的停止准则直接相关。

2.2 随机欠采样如何构造测量矩阵

压缩感知的第二个前提是非相干测量，也就是测量矩阵 Φ（大小 M×N，M 远小于 N）与稀疏基 Ψ（此处为 DFT 矩阵）之间要满足受限等距性质。常见做法是直接用随机高斯矩阵或抽取部分行实现欠采样。对正弦信号最直观的欠采样方式，是保留原始等间隔采样信号的少量随机位置样本，这相当于 Φ 是稀疏行抽取矩阵，实现成本最低。

rng(42); % 固定随机种子,保证可复现 M = 120; % 欠采样点数,约为 N/4 idx = sort(randperm(N, M)); % 从 512 个点中随机抽取 120 个位置 y = x(idx); % 欠采样后的观测值

代码说明：randperm(N, M)生成 M 个不重复的随机索引，sort保证时间顺序排列，观测向量 y 就是这些随机时间点上的信号幅值。随机种子的设置决定了每次运行的测量位置，建议在对比实验时固定同一个种子，避免测量矩阵不同导致结论失真。

测量矩阵的另一种更标准的写法是直接构造高斯随机矩阵 Φ，与稀疏基 Ψ 相乘得到感知矩阵 A，这在恢复时用得更普遍。

Phi = randn(M, N)/sqrt(M); % 高斯随机测量矩阵 Psi = dftmtx(N)/sqrt(N); % 归一化 DFT 稀疏基 A = Phi * Psi; % 感知矩阵 M×N y = Phi * x(:); % 观测向量

参数说明：高斯矩阵乘以系数1/sqrt(M)是为了让 Φ 的行范数接近 1，避免恢復算法的阈值设置依赖信号尺度；Psi是归一化 DFT 基，保证变换系数幅值与信号幅度同一量级，后续 OMP 的残差阈值也好设定。

2.3 恢复算法选型：凸优化与贪婪算法的取舍

有了观测向量 y 和感知矩阵 A，问题变成求解欠定方程 y = A·s，其中 s 是信号在 DFT 域的系数向量。由于 K 远小于 N，求解可以转化为最小化 l0 范数的组合优化问题，但 l0 是 NP 困难问题，工程上走两条路：基追踪（Basis Pursuit）把它松弛为 l1 凸优化，用 CVX 或 SPGL1 求解；另一种是系列贪婪算法，以正交匹配追踪（OMP）为代表，每次迭代选一个与残差最相关的原子，逐步逼近真实支撑集。

我一般在 MATLAB 里优先用 OMP 做教学和快速验证，因为它实现简单、迭代过程透明，稀疏度 K 已知时效果好；基追踪的优势在于不精确已知 K 也能用，适合噪声环境，但安装 CVX 对新手门槛高。两种算法在正弦恢复上的效果差距通常在 1dB 以内，所以从 OMP 入门是效率最高的路径。

3. MATLAB 实现多正弦信号随机欠采样的完整流程与参数选择

3.1 构造多频正弦信号的注意事项

实际仿真中叠加正弦的频率不能随意选，要保证在 DFT 域确实是稀疏的，也就是频率落在 DFT 频点上或接近频点。如果选择非整周期频率如 53.7Hz，FFT 会泄露到附近很多频点，稀疏度迅速上升，压缩感知的前提就不成立，恢复误差大幅增加。

N = 512; fs = 1024; t = (0:N-1)/fs; % 总时长 0.5s f_set = [50, 128, 257]; % 注意 128 和 257 都是整周期频率 amp_set = [1.0, 0.6, 0.3]; phase_set = [0, pi/4, pi/3]; x = zeros(1, N); for k = 1:3 x = x + amp_set(k) * sin(2*pi*f_set(k)*t + phase_set(k)); end

说明：频率设为整周期是为了避免频谱泄露。128 对应 64 个完整周期，257 接近奈奎斯特频域但仍在范围内，此时 DFT 的峰值恰好落在对应的频点上，非零系数只有 3 个（K=3）。相位随机化不影响稀疏度，但能验证算法对相位不敏感的性质。

3.2 随机欠采样点的个数与恢复质量的关系

欠采样点数 M 直接决定测量数，理论上 M ≥ C·K·log(N/K) 就能大概率精确恢复，C 是常数，经验上取 2～4。对 N=512、K=3 的情况，理论下限大约 20 个点需要，但考虑到数值稳定性，M 取 80~150 是合理区间。取点的位置如果是纯随机，可能出现局部间隔过大的极端分布，对恢复带来不确定性，所以可以采用分段随机的方式，比如把整个时间轴平分成若干段，每段内随机抽一个或两个点。

seg = 8; % 等分为 8 段 pts_per_seg = 15; % 每段取 15 个点,共 120 点 idx = []; for k = 1:seg seg_start = (k-1)*N/seg + 1; seg_end = k*N/seg; idx_seg = sort(randperm(N/seg, pts_per_seg)) + seg_start - 1; idx = [idx, idx_seg]; end y = x(idx);

这种分段随机策略比全局随机更稳定，避免极端情况下全部采样点挤在前半段导致后半段信息割裂。如果采用高斯随机测量矩阵，则不是抽取时间点而是对全部信号做随机投影，不依赖采样位置的均匀性，但计算量大一些。

3.3 欠采样前是否需要加抗混叠滤波器

工程上直接对连续信号做随机欠采样时，模拟端通常要加带宽限制滤波器，否则高频成分会混叠到低频造成信号污染。但在仿真层面，信号本来就是在数字域构造的，所以不存在真实的混叠过程，只需要在构造信号时确保最高频率低于 f_s/2。需要注意的是，如果随机欠采样后的等效采样率低于奈奎斯特率，信号在某些时间段看似丢失了高频信息，这正是压缩感知要解决的问题——利用频域稀疏性把这些信息重新补回来。

4. OMP 恢复算法在 MATLAB 中的实现与参数调试

4.1 正交匹配追踪的迭代逻辑与 MATLAB 代码

OMP 的核心逻辑分四步：计算感知矩阵各列与残差的相关系数，选出最相关的一列加入支撑集，用最小二乘估计系数，更新残差并继续迭代，直到达到稀疏度或残差阈值。每步都要保证支撑集列之间尽量正交，所以叫正交匹配追踪。

function s_hat = my_omp(A, y, K) [M, N] = size(A); r = y; % 残差初始化 idx_selected = zeros(K, 1); % 记录被选列序号 A_selected = zeros(M, K); % 记录被选列 s_hat = zeros(N, 1); for k = 1:K % 1. 计算残差与每个原子的相关系数 corr = abs(A' * r); % 2. 选择相关系数最大的原子下标 [~, max_idx] = max(corr); idx_selected(k) = max_idx; A_selected(:, k) = A(:, max_idx); % 3. 最小二乘求解当前支撑集下的系数 A_sub = A_selected(:, 1:k); coef = A_sub \ y; % 4. 更新残差 r = y - A_sub * coef; end s_hat(idx_selected) = coef; end

这段代码有几处关键点需要说明。每次迭代更新系数时是对所有已选原子一起做最小二乘，而不是只更新最新原子的系数，因为这个步骤保证了残差始终与已选原子张成的空间正交。对于 K=3，每轮A_sub \ y的求解规模不超过 512×3，极快，即使 N 增大到上万也能接受。需要留意max(corr)遇到并列时会选第一个位置，如果信号本身有两个频率完全对称的分量，并列概率会上升，此时可以给相关性加微小的随机扰动打破并列。

4.2 从恢复系数反变换回时域信号

OMP 解出的是频域系数向量 s_hat，要得到恢复信号还需要左乘稀疏基的逆矩阵。由于稀疏基是归一化 DFT 矩阵，它的逆就是本身的共轭转置。

x_rec = Psi' * s_hat; % x_rec 是时域恢复信号 freq_axis = (0:N-1)/N*fs; figure; subplot(2,1,1); plot(t, x, 'b', t, x_rec, 'r--'); subplot(2,1,2); stem(abs(s_hat), 'filled');

接着计算相对误差，公式为err = norm(x(:)-x_rec(:)) / norm(x(:))。误差低于 1e-6 说明恢复完全正确；误差在 1e-2 量级说明系数匹配有问题，大概率是频率设定落在了非整周期频点上；误差在 0.1 以上则基本是算法参数或测量点数的问题。频谱图上如果恢复结果的谱线出现在正确位置但幅度偏差，优先检查最小二乘步骤有没有把幅值正确映射回，而不是怀疑算法本身。

4.3 OMP 参数调节的优先级顺序

参数调节的最有效路径是：先固定稀疏度K，调节测量点数M，看恢复误差的拐点；然后固定M，把K从1到6变化观察误差变化规律，如果K设得比真实值小，误差会居高不下，如果K设得偏大，多余迭代只在数值噪声里打转，误差略升但不会灾难性恶化；最后调噪声环境下的停止阈值。

Klist = 1:6; errlist = zeros(size(Klist)); for i = 1:length(Klist) s_hat = my_omp(A, y, Klist(i)); x_rec = Psi' * s_hat; errlist(i) = norm(x(:)-x_rec(:)) / norm(x(:)); end plot(Klist, errlist);

这段代码把稀疏度从1遍历到6，观察误差曲线在K=3处是否有显著拐点。实践中这是最常用的调试手段，也是验证算法正确性的第一步。因为真实稀疏度已知，曲线的形态应该足够明显——K从1到2和从2到3误差逐步下降，跨过3之后误差基本维持在 1e-6 水平。如果曲线在K=3处没有明显拐点，说明感知矩阵的计算过程有错误，优先检查 A = Phi*Psi 的构造是否少乘了系数。

5. 测量矩阵的进阶对比与噪声鲁棒性分析

5.1 高斯矩阵、伯努利矩阵与随机抽取矩阵的实测对比

上面用的是随机抽取部分时间样本的方式，接下来直接比较三种常见测量矩阵：高斯随机矩阵、伯努利±1矩阵、随机抽取单位阵行。为公平对比，保持M=120、K=3、N=512不变，用同一信号分别做恢复。

% 高斯矩阵 Phi_gauss = randn(M, N)/sqrt(M); % 伯努利矩阵 Phi_bern = (rand(M, N) > 0.5)*2 - 1; % ±1 Phi_bern = Phi_bern / sqrt(M); % 抽取矩阵 idx_rand = randperm(N, M); Phi_sample = zeros(M, N); for i = 1:M Phi_sample(i, idx_rand(i)) = 1; end

随机抽取矩阵本质上就是时间欠采样，恢复误差由采样位置的分布决定。高斯矩阵对任意稀疏信号都有良好的非相干性，但计算量最大，因为它要和DFT矩阵做乘法生成 A，M×N 乘法的开销在数据规模增大时线性上涨。伯努利矩阵的优势是存储和计算更快，但信号本身如果是极稀疏的，效果和高斯相当。实测中，三种矩阵在K=3时恢复误差都能到1e-6量级，真正的区别体现在高稀疏度或强噪声场景下——高斯矩阵的稳定性最好，伯努利次之，抽取矩阵最不稳定。

5.2 恢复误差随采样率变化的模拟

采样率从 N/8 逐步涨到 N/2，记录恢复误差的对数值。

M_ratio = [0.1, 0.15, 0.2, 0.25, 0.3, 0.4, 0.5]; errs = zeros(size(M_ratio)); for i = 1:length(M_ratio) M_cur = round(N * M_ratio(i)); Phi_cur = randn(M_cur, N)/sqrt(M_cur); A_cur = Phi_cur * Psi; y_cur = Phi_cur * x(:); s_hat = my_omp(A_cur, y_cur, 3); x_rec = Psi' * s_hat; errs(i) = norm(x(:)-x_rec(:)) / norm(x(:)); end semilogy(M_ratio, errs, 'o-');

运行结果通常呈现一条陡峭下降的曲线：采样率低于某个阈值时误差在 0.1 以上，超过阈值后迅速跌到 1e-6 以下。这个阈值附近就是相变区，代表最低所需采样率。工程估值不要追求在相变区内工作，因为数值条件差，微小的噪声或参数扰动都会让恢复结果大幅波动，设计时应把采样率取在相变区右侧至少 1.2 倍的位置。这里的相变点直接由 K/N 决定，K 越大相变点越靠右。

5.3 有噪环境下的参数调整

真实场景里观测信号总是带噪声的，此时 y = Phi*x + n，n 是高斯白噪声。OMP 在噪声下容易出现两个问题：一是残差阈值设得太严，把噪声当信号继续迭代；二是最小二乘阶段对噪声敏感，系数估计方差增大。针对第一个问题，可以改用相对残差下降作为停止条件，即当norm(r) < 1e-6 * norm(y)时停止迭代，避免迭代次数过多拟合噪声。针对第二个问题，可以在最小二乘步骤加一个小的 Tikhonov 正则项。

% 带正则的最小二乘 lambda = 1e-3; coef = (A_sub'*A_sub + lambda*eye(k)) \ (A_sub'*y);

正则系数 lambda 的选择有个经验范围：信噪比 20dB 以下时取 1e-2，信噪比 40dB 以上时取 1e-4。lambda 太大把真实系数也压平了，太小起不到正则作用。以 20dB 噪声为例，不加速正时恢复误差约 0.08，加正则后可以降到 0.03 左右，效果直观可见。另外，噪声环境下建议改用基追踪或 SPGL1 这类凸优化算法，它们在噪声处理上理论保证更强，MATLAB 里可以调用 SPGL1 工具箱函数，内部实现了谱投影梯度求解。

6. 一小时内复现完整压缩感知正弦信号恢复的验证技巧

6.1 组合全部步骤的可运行脚本

把以上各段拼成一个完整脚本cs_sin_demo.m，从生成信号到展示恢复结果的完整代码如下，直接保存运行即可看到恢复效果。

%% 生成多正弦信号 N = 512; fs = 1024; t = (0:N-1)/fs; f_set = [50, 128, 257]; x = sin(2*pi*f_set(1)*t) + 0.6*sin(2*pi*f_set(2)*t + pi/4) + 0.3*sin(2*pi*f_set(3)*t + pi/3); %% 随机欠采样测量 rng(42); M = 120; Phi = randn(M, N)/sqrt(M); y = Phi * x(:); %% 稀疏基和感知矩阵 Psi = dftmtx(N)/sqrt(N); A = Phi * Psi; %% OMP 恢复 s_hat = my_omp(A, y, 3); x_rec = Psi' * s_hat; %% 误差与可视化 err = norm(x(:)-x_rec(:)) / norm(x(:)); fprintf('恢复相对误差: %.2e\n', err); figure; subplot(2,1,1); plot(t, x, 'b'); hold on; plot(t, real(x_rec), 'r--'); legend('原始', '恢复'); subplot(2,1,2); stem(abs(s_hat)); title('恢复的频域系数');

运行这个脚本后，第一行输出应该是恢复相对误差: 1e-15左右的数值，如果电脑精度较高可能显示为 0。这个结果是验证算法实现正确的最快途径。如果误差很大，逐个环节排查：先检查 y 的数值范围，再用cond(A'*A)检查感知矩阵的条件数。

6.2 用一致性校验判断恢复是否可信

工程上有个不依赖已知信号的校验技巧：把恢复信号重新作为原始信号再做一次随机欠采样，第二次的观测值和第一次对不上，说明恢复结果与真实信号不一致，算法过程有问题或测量点数不足。这个技巧在无法预先知道信号真实值时特别有用，它本质上是把压缩感知问题当作一个可重入的编解码回路来验证。

Phi2 = randn(M, N)/sqrt(M); y2 = Phi2 * x_rec(:); y2_true = Phi2 * real(x_rec(:)); consistency = norm(y2 - y2_true) / norm(y2);

当一致性偏差远大于恢复误差时，优先怀疑 OMP 的迭代停止条件没有满足，或感知矩阵 A 的构造步骤出现错误。反之如果一致性不比恢复误差高太多，说明恢复结果基本可信。

6.3 压缩感知恢复正弦信号的三处易错点

第一处是忘记归一化 DFT 矩阵，直接用未归一化的dftmtx(N)，这会让稀疏系数尺度偏离真实值，稀疏度判断出错。第二处是 OMP 的最小二乘没有使用全支撑集的共同求解，而只是更新最新一个原子的系数，这一步导致残差更新不彻底。第三处是随机测量矩阵没有固定随机种子，同一段代码前后跑两次结果不同，无法复现实验结论。检查完这三处，剩余的调试路径就是稀疏度和测量点数了。

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