简介一份以压缩感知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 三个位置凸起其余位置接近机器精度。K3 就是信号的稀疏度也是后面测量矩阵设计、恢复算法迭代次数设定的核心依据。稀疏度估计偏低时恢复误差迅速恶化偏高时计算量浪费但结果仍然正确这与匹配追踪类算法的停止准则直接相关。2.2 随机欠采样如何构造测量矩阵压缩感知的第二个前提是非相干测量也就是测量矩阵 Φ大小 M×NM 远小于 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.7HzFFT 会泄露到附近很多频点稀疏度迅速上升压缩感知的前提就不成立恢复误差大幅增加。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 个K3。相位随机化不影响稀疏度但能验证算法对相位不敏感的性质。3.2 随机欠采样点的个数与恢复质量的关系欠采样点数 M 直接决定测量数理论上 M ≥ C·K·log(N/K) 就能大概率精确恢复C 是常数经验上取 24。对 N512、K3 的情况理论下限大约 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这段代码有几处关键点需要说明。每次迭代更新系数时是对所有已选原子一起做最小二乘而不是只更新最新原子的系数因为这个步骤保证了残差始终与已选原子张成的空间正交。对于 K3每轮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观察误差曲线在K3处是否有显著拐点。实践中这是最常用的调试手段也是验证算法正确性的第一步。因为真实稀疏度已知曲线的形态应该足够明显——K从1到2和从2到3误差逐步下降跨过3之后误差基本维持在 1e-6 水平。如果曲线在K3处没有明显拐点说明感知矩阵的计算过程有错误优先检查 A Phi*Psi 的构造是否少乘了系数。5. 测量矩阵的进阶对比与噪声鲁棒性分析5.1 高斯矩阵、伯努利矩阵与随机抽取矩阵的实测对比上面用的是随机抽取部分时间样本的方式接下来直接比较三种常见测量矩阵高斯随机矩阵、伯努利±1矩阵、随机抽取单位阵行。为公平对比保持M120、K3、N512不变用同一信号分别做恢复。% 高斯矩阵 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矩阵做乘法生成 AM×N 乘法的开销在数据规模增大时线性上涨。伯努利矩阵的优势是存储和计算更快但信号本身如果是极稀疏的效果和高斯相当。实测中三种矩阵在K3时恢复误差都能到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 nn 是高斯白噪声。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 的最小二乘没有使用全支撑集的共同求解而只是更新最新一个原子的系数这一步导致残差更新不彻底。第三处是随机测量矩阵没有固定随机种子同一段代码前后跑两次结果不同无法复现实验结论。检查完这三处剩余的调试路径就是稀疏度和测量点数了。本文还有配套的精品资源点击获取