在信号处理领域，我们经常面对这样的问题：传感器采集到的信号往往含有大量噪声，而信号本身可能是稀疏的（即大部分元素为零或接近零）。如何从含噪观测中精确恢复原始稀疏信号？这不仅是压缩感知的核心问题，也是许多工程应用（如医学成像、雷达探测、地震勘探）的关键。什么是稀疏信号？
稀疏信号 是指信号中只有少数元素非零，其余都为零（或可忽略不计）。例如：
一个由几个脉冲组成的信号
图像经过小波变换后，大部分系数接近于零
语音信号在特定字典下的表示系数
数学上，如果一个长度为 N 的信号 x只有 K 个非零元素（K《N），我们就称它是 K-稀疏的。
但在实际中，我们观测到的信号往往被噪声污染，比如传感器热噪声、量化噪声等。此时，我们得到的是含噪信号 y=x+n，其中 n 为噪声。如何从 y 中尽可能恢复出干净的 x？这便是稀疏信号去噪的任务。
稀疏恢复的核心思想：利用稀疏性作为先验
如果对信号没有任何先验知识，从含噪观测中恢复原信号是病态的。但如果我们知道信号是稀疏的，就可以利用这一先验来约束解空间。
1. 基于L1范数的优化方法
最经典的稀疏恢复模型是基追踪去噪（Basis Pursuit Denoising，BPDN）：
2. 贪婪算法：正交匹配追踪（OMP）
另一种直观的方法是正交匹配追踪（Orthogonal Matching Pursuit，OMP）。它迭代地选择与当前残差最相关的原子，然后通过最小二乘更新系数，直到满足停止条件（如稀疏度或残差阈值）。OMP简单高效，非常适合理解稀疏恢复的流程。
Python实战：用OMP从噪声中恢复稀疏信号
import numpy as npimport matplotlib.pyplot as pltimport seaborn as sns# 设置绘图风格sns.set_style("whitegrid")plt.rcParams['font.sans-serif'] = ['SimHei'] # 显示中文plt.rcParams['axes.unicode_minus'] = Falsedef generate_sparse_signal(n, k): """生成一个长度为n，具有k个非零值的稀疏信号""" x = np.zeros(n) pos = np.random.choice(n, k, replace=False) # 随机选择非零位置 x[pos] = np.random.randn(k) * 3 # 随机幅度 return xdef add_noise(x, snr_db): """添加高斯白噪声，SNR单位为dB""" signal_power = np.mean(x**2) noise_power = signal_power / (10**(snr_db/10)) noise = np.sqrt(noise_power) * np.random.randn(len(x)) return x + noisedef omp(y, k, max_iter=None): """ 正交匹配追踪算法 参数： y: 观测信号（含噪） k: 稀疏度（非零元素个数） max_iter: 最大迭代次数（默认等于k） 返回： x_hat: 恢复的稀疏信号 """ n = len(y) if max_iter is None: max_iter = k x_hat = np.zeros(n) residual = y.copy() support = [] # 支撑集（已选原子索引） for _ in range(max_iter): # 计算与残差的相关性 correlations = np.abs(residual) # 这里我们使用标准基（单位矩阵），所以相关就是残差本身 # 实际上对于一般字典 D，应为 D.T @ residual，此处D=I # 选择相关性最大的索引（注意已选的不再选） j = np.argmax(correlations) if j in support: # 避免重复选择（如果已经选过，选择次大的） correlations[j] = -np.inf j = np.argmax(correlations) support.append(j) # 用当前支撑集通过最小二乘求解系数 # 由于字典是标准基，支撑集上的系数就是对应位置的y值 # 但为了通用性，我们显式构造子字典 D_s = I[:, support] D_s = np.eye(n)[:, support] # 子字典，大小为 n x |support| # 最小二乘解：x_s = (D_s^T D_s)^{-1} D_s^T y # 因为 D_s 的列是标准正交的（单位矩阵的列），所以 D_s^T D_s = I，因此 x_s = D_s^T y x_s = D_s.T @ y # 实际上就是 y[support] # 更新残差 residual = y - D_s @ x_s # 可选：检查残差是否足够小 if np.linalg.norm(residual) < 1e-6: break # 将恢复的系数放回全长向量 x_hat[support] = x_s return x_hat# ================== 参数设置 ==================n = 256 # 信号长度k = 15 # 稀疏度（非零个数）snr_db = 10 # 信噪比 (dB)# ================== 生成数据 ==================np.random.seed(42) # 固定随机种子，便于复现x_true = generate_sparse_signal(n, k)y_noisy = add_noise(x_true, snr_db)# ================== OMP恢复 ==================x_omp = omp(y_noisy, k)# ================== 结果可视化 ==================fig, axes = plt.subplots(3, 1, figsize=(12, 8), sharex=True)# 原始稀疏信号axes[0].stem(x_true, linefmt='C0-', markerfmt='C0o', basefmt='k-', use_line_collection=True)axes[0].set_ylabel('幅度', fontsize=12)axes[0].set_title('原始稀疏信号 (k={})'.format(k), fontsize=14)# 含噪信号axes[1].plot(y_noisy, color='C1', linewidth=1)axes[1].set_ylabel('幅度', fontsize=12)axes[1].set_title('含噪观测 (SNR={} dB)'.format(snr_db), fontsize=14)# OMP恢复信号axes[2].stem(x_omp, linefmt='C2-', markerfmt='C2o', basefmt='k-', use_line_collection=True)axes[2].set_xlabel('样本点', fontsize=12)axes[2].set_ylabel('幅度', fontsize=12)axes[2].set_title('OMP恢复信号 (稀疏度 k={})'.format(k), fontsize=14)plt.tight_layout()plt.show()# 计算恢复误差mse = np.mean((x_true - x_omp)**2)print(f"均方误差 (MSE): {mse:.6f}")
