发布于 由admin

Python实现基于压缩感知（Compressed Sensing）的稀疏采样与信号恢复

Python实现基于压缩感知（Compressed Sensing）的稀疏采样与信号恢复

大家好，今天我们来探讨一个在信号处理领域非常重要的技术——压缩感知（Compressed Sensing，CS）。CS的核心思想是在远低于奈奎斯特采样率的情况下，对稀疏或可压缩信号进行采样，并利用特定的算法精确或近似地恢复原始信号。这在数据采集、图像处理、医学成像等领域都有着广泛的应用前景。

本次讲座将以Python编程为基础，深入讲解CS的原理，并提供相应的代码实现，帮助大家理解并掌握这一技术。

1. 压缩感知的基本原理

传统的奈奎斯特-香农采样定理指出，为了无失真地恢复信号，采样频率至少要大于信号最高频率的两倍。然而，在很多实际应用中，信号本身是稀疏的，即信号的大部分能量集中在少数几个频率分量上。压缩感知正是利用了这一特性。

1.1 稀疏性

一个信号 x ∈ ℝN 被称为 K-稀疏的，如果它最多只有 K 个非零元素，其中 K << N。在实际应用中，信号可能不是严格稀疏的，但可以在某个变换域（例如，傅里叶变换、小波变换）下变得稀疏或近似稀疏。

1.2 测量矩阵

压缩感知不是直接采样原始信号 x，而是通过一个 M × N 的测量矩阵 Φ 对信号进行线性测量，得到测量值 y ∈ ℝM，其中 M < N：

y = Φx

关键在于，即使 M << N，我们仍然可以通过合适的算法从 y 中恢复出原始信号 x。

1.3 恢复算法

压缩感知的核心在于找到一个有效的恢复算法，从测量值 y 和测量矩阵 Φ 中重构出稀疏信号 x。常见的恢复算法包括：

  • 凸优化方法： 将稀疏信号的恢复问题转化为一个凸优化问题，例如，L1范数最小化。
  • 贪婪算法： 通过迭代的方式，逐步逼近稀疏信号，例如，匹配追踪（Matching Pursuit，MP）、正交匹配追踪（Orthogonal Matching Pursuit，OMP）。

2. Python实现：生成稀疏信号与测量

首先，我们需要生成一个稀疏信号和一个测量矩阵。

import numpy as np
import matplotlib.pyplot as plt

def generate_sparse_signal(N, K):
  """生成一个长度为N，稀疏度为K的稀疏信号"""
  x = np.zeros(N)
  indices = np.random.choice(N, K, replace=False)  # 随机选择K个位置
  x[indices] = np.random.randn(K)  # 在这些位置上赋予随机值
  return x

def generate_measurement_matrix(M, N):
  """生成一个M x N的测量矩阵，使用高斯随机矩阵"""
  Phi = np.random.randn(M, N)
  return Phi

# 设置参数
N = 256  # 信号长度
K = 32   # 稀疏度
M = 64   # 测量次数

# 生成稀疏信号和测量矩阵
x = generate_sparse_signal(N, K)
Phi = generate_measurement_matrix(M, N)

# 进行测量
y = np.dot(Phi, x)

# 可视化
plt.figure(figsize=(12, 6))

plt.subplot(1, 2, 1)
plt.stem(x)
plt.title("Sparse Signal")
plt.xlabel("Index")
plt.ylabel("Amplitude")

plt.subplot(1, 2, 2)
plt.stem(y)
plt.title("Measurements")
plt.xlabel("Index")
plt.ylabel("Amplitude")

plt.tight_layout()
plt.show()

这段代码首先定义了两个函数：generate_sparse_signal 用于生成一个 K-稀疏的信号，generate_measurement_matrix 用于生成一个高斯随机测量矩阵。然后，我们设置信号长度 N、稀疏度 K 和测量次数 M，并使用这两个函数生成稀疏信号 x 和测量矩阵 Φ。最后，我们通过 y = Φx 得到测量值 y，并使用Matplotlib将原始稀疏信号和测量值进行可视化。

3. Python实现：基于L1范数最小化的信号恢复

接下来，我们将使用凸优化方法，通过L1范数最小化来恢复稀疏信号。这里我们使用 cvxopt 库来解决凸优化问题。

from cvxopt import matrix, solvers

def l1_recovery(Phi, y, N):
  """使用L1范数最小化恢复稀疏信号"""
  M = Phi.shape[0]

  # 定义优化问题的系数矩阵
  c = matrix(np.concatenate((np.zeros(N), np.ones(M)), axis=0))
  G = matrix(np.vstack((np.hstack((Phi, -np.eye(M))),
                        np.hstack((-Phi, -np.eye(M))))))
  h = matrix(np.concatenate((y, -y), axis=0))

  # 求解线性规划问题
  solvers.options['show_progress'] = False  # 静默模式，不显示solver的迭代过程
  solution = solvers.lp(c, G, h)

  # 提取恢复的信号
  x_recovered = np.array(solution['x'][:N]).flatten()
  return x_recovered

# 恢复信号
x_recovered = l1_recovery(Phi, y, N)

# 可视化
plt.figure(figsize=(12, 6))

plt.subplot(1, 2, 1)
plt.stem(x)
plt.title("Original Sparse Signal")
plt.xlabel("Index")
plt.ylabel("Amplitude")

plt.subplot(1, 2, 2)
plt.stem(x_recovered)
plt.title("Recovered Signal (L1 Minimization)")
plt.xlabel("Index")
plt.ylabel("Amplitude")

plt.tight_layout()
plt.show()

# 评估恢复效果
error = np.linalg.norm(x - x_recovered) / np.linalg.norm(x)
print(f"Relative Error: {error}")

这段代码定义了一个 l1_recovery 函数，该函数使用 cvxopt 库来解决L1范数最小化问题。该函数将压缩感知的恢复问题转化为一个线性规划问题，并使用 cvxopt.solvers.lp 函数进行求解。最后，提取恢复的信号并返回。为了方便观察效果，我们计算了相对误差。

解释：

L1范数最小化问题的数学表达形式为：

minimize ||x||1 subject to Φx = y

其中 ||x||1 表示 x 的L1范数，即所有元素的绝对值之和。

在实际求解时，我们通常将其转化为一个等价的线性规划问题：

minimize ∑ ti subject to –ti ≤ xi ≤ ti and Φx = y

为了能被 cvxopt 库处理，我们将此问题转化为标准形式：

minimize cTz subject to Gz ≤ h

其中 z = [x, t]T，c, G, h 需要根据原问题进行构造。

4. Python实现：基于正交匹配追踪（OMP）的信号恢复

除了L1范数最小化，我们还可以使用贪婪算法，例如正交匹配追踪（OMP）来恢复稀疏信号。

def orthogonal_matching_pursuit(Phi, y, K):
  """使用正交匹配追踪算法恢复稀疏信号"""
  N = Phi.shape[1]
  residual = y.copy()
  indices = []
  x_recovered = np.zeros(N)

  for _ in range(K):
    # 找到与残差相关性最大的原子
    correlations = np.abs(np.dot(Phi.T, residual))
    index = np.argmax(correlations)
    indices.append(index)

    # 使用最小二乘法更新信号
    Phi_subset = Phi[:, indices]
    x_subset = np.linalg.lstsq(Phi_subset, y, rcond=None)[0]

    # 更新恢复的信号和残差
    x_recovered[indices] = x_subset
    residual = y - np.dot(Phi, x_recovered)

  return x_recovered

# 恢复信号
x_recovered_omp = orthogonal_matching_pursuit(Phi, y, K)

# 可视化
plt.figure(figsize=(12, 6))

plt.subplot(1, 2, 1)
plt.stem(x)
plt.title("Original Sparse Signal")
plt.xlabel("Index")
plt.ylabel("Amplitude")

plt.subplot(1, 2, 2)
plt.stem(x_recovered_omp)
plt.title("Recovered Signal (OMP)")
plt.xlabel("Index")
plt.ylabel("Amplitude")

plt.tight_layout()
plt.show()

# 评估恢复效果
error_omp = np.linalg.norm(x - x_recovered_omp) / np.linalg.norm(x)
print(f"Relative Error (OMP): {error_omp}")

这段代码定义了一个 orthogonal_matching_pursuit 函数，该函数实现了 OMP 算法。OMP 算法通过迭代的方式，每次选择与残差相关性最大的原子（测量矩阵的列），然后使用最小二乘法更新信号，并更新残差。该算法重复执行 K 次，其中 K 是信号的稀疏度。

OMP算法步骤：

  1. 初始化： 残差 r = y, 索引集合 Λ = ∅, 恢复的信号 x = 0。
  2. 迭代：
    • 找到与残差相关性最大的原子：λ = argmaxj |<r, Φj>|
    • 更新索引集合：Λ = Λ ∪ {λ}
    • 使用最小二乘法更新信号：xΛ = (ΦΛTΦΛ)-1 ΦΛTy
    • 更新残差：r = y – ΦΛxΛ
  3. 停止： 当迭代次数达到 K 时，停止迭代。

5. 不同恢复算法的比较

为了更直观地比较 L1 范数最小化和 OMP 算法的性能，我们可以使用不同的测量次数 M，并计算恢复误差。

M_values = [32, 64, 128]  # 不同的测量次数
errors_l1 = []
errors_omp = []

for M in M_values:
  # 生成测量矩阵
  Phi = generate_measurement_matrix(M, N)
  y = np.dot(Phi, x)

  # 使用L1范数最小化恢复
  x_recovered_l1 = l1_recovery(Phi, y, N)
  error_l1 = np.linalg.norm(x - x_recovered_l1) / np.linalg.norm(x)
  errors_l1.append(error_l1)

  # 使用OMP恢复
  x_recovered_omp = orthogonal_matching_pursuit(Phi, y, K)
  error_omp = np.linalg.norm(x - x_recovered_omp) / np.linalg.norm(x)
  errors_omp.append(error_omp)

# 绘制误差曲线
plt.figure(figsize=(8, 6))
plt.plot(M_values, errors_l1, marker='o', label='L1 Minimization')
plt.plot(M_values, errors_omp, marker='o', label='OMP')
plt.xlabel("Number of Measurements (M)")
plt.ylabel("Relative Error")
plt.title("Comparison of Recovery Algorithms")
plt.legend()
plt.grid(True)
plt.show()

# 将结果以表格形式展示
data = {'M': M_values,
        'L1 Error': errors_l1,
        'OMP Error': errors_omp}

print("nRecovery Error Comparison:")
print(pd.DataFrame(data))

这段代码首先定义了不同的测量次数 M 的值，然后针对每个 M 值，生成测量矩阵并进行测量，分别使用 L1 范数最小化和 OMP 算法进行信号恢复，并计算恢复误差。最后，将结果以图表和表格的形式进行展示，方便我们比较不同算法的性能。

6. 讨论：影响压缩感知性能的因素

压缩感知的性能受到多种因素的影响，包括：

  • 信号的稀疏度 K： 信号越稀疏，越容易被恢复。
  • 测量矩阵 Φ 的选择： 测量矩阵需要满足RIP（Restricted Isometry Property）性质，以保证信号可以被精确或近似恢复。高斯随机矩阵是一种常用的测量矩阵。
  • 测量次数 M： 测量次数越多，恢复效果越好，但同时也会增加数据采集的成本。
  • 恢复算法的选择： 不同的恢复算法具有不同的性能和计算复杂度。
  • 噪声： 噪声会降低恢复效果，需要使用鲁棒的恢复算法。

7. 更深入的优化和应用

上面的代码仅仅是演示压缩感知基本原理的简单例子。在实际应用中，需要根据具体情况进行优化，例如：

  • 选择更合适的测量矩阵： 例如 Toeplitz 矩阵、循环矩阵等，可以降低存储和计算复杂度。
  • 使用更高级的恢复算法： 例如，迭代硬阈值算法（Iterative Hard Thresholding，IHT）、 approximate message passing (AMP) 等，可以提高恢复速度和精度。
  • 考虑信号的先验信息： 例如，信号的非负性、有界性等，可以提高恢复效果。
  • 在不同领域中的应用： 例如，图像处理、医学成像、无线通信等，需要根据具体应用场景进行定制。
import pywt  # 导入小波变换库

def wavelet_transform(signal, wavelet='db4', level=3):
    """
    对信号进行小波变换.
    """
    coeffs = pywt.wavedec(signal, wavelet, level=level)
    return coeffs

def inverse_wavelet_transform(coeffs, wavelet='db4'):
    """
    从小波系数重构信号.
    """
    signal = pywt.waverec(coeffs, wavelet)
    return signal

def sparsify_wavelet(signal, K):
    """
    小波域稀疏化：保留最大的 K 个系数，其余置零.
    """
    coeffs = wavelet_transform(signal)
    all_coeffs = np.concatenate([c.ravel() for c in coeffs])
    indices = np.argsort(np.abs(all_coeffs))[::-1]  # 绝对值降序排列的索引
    mask = np.zeros_like(all_coeffs, dtype=bool)
    mask[indices[:K]] = True  # 保留最大的 K 个系数

    start = 0
    sparse_coeffs = []
    for c in coeffs:
        end = start + c.size
        sparse_coeffs.append(all_coeffs[start:end] * mask[start:end].astype(float).reshape(c.shape))
        start = end

    return sparse_coeffs

#示例使用：对信号进行小波变换，保留最大的K个系数，然后重构信号
# 生成一个测试信号
N = 256
signal = np.random.randn(N)

# 小波变换并稀疏化
K = 64 # 保留64个最大的系数
sparse_coeffs = sparsify_wavelet(signal, K)

# 重构信号
reconstructed_signal = inverse_wavelet_transform(sparse_coeffs)

# 可视化结果
plt.figure(figsize=(12, 6))

plt.subplot(1, 2, 1)
plt.plot(signal)
plt.title("Original Signal")
plt.xlabel("Index")
plt.ylabel("Amplitude")

plt.subplot(1, 2, 2)
plt.plot(reconstructed_signal)
plt.title("Reconstructed Signal (Wavelet Sparsity)")
plt.xlabel("Index")
plt.ylabel("Amplitude")

plt.tight_layout()
plt.show()

这个示例演示了如何利用小波变换来提高压缩感知的性能。首先，对信号进行小波变换，然后保留最大的 K 个小波系数，将其他系数置零，从而实现信号的稀疏化。最后，使用逆小波变换重构信号。

结论： 压缩感知为信号处理提供了新思路

本次讲座我们深入探讨了压缩感知的基本原理和 Python 实现。从生成稀疏信号到使用 L1 范数最小化和 OMP 算法进行信号恢复，再到不同恢复算法的比较，以及影响压缩感知性能的因素的讨论，相信大家对压缩感知有了更清晰的理解。

压缩感知不仅是一种信号处理技术，更是一种新的数据采集和处理的思路。它在数据量爆炸式增长的今天，具有重要的应用价值。希望通过本次讲座，能够激发大家对压缩感知的兴趣，并在实际应用中发挥它的潜力。

更多IT精英技术系列讲座，到智猿学院

发表回复

您的邮箱地址不会被公开。 必填项已用 * 标注