压缩感知入门与 Matlab 教程

引言

压缩感知 (CS) 是信号处理领域一项相对较新的技术，它允许在仅采集少量样本的情况下获取信号。它适用于稀疏信号，并有一些限制，我们将在后续介绍。

对于熟悉奈奎斯特定理的人来说，它指出为了获取信号中的所有相关信息，采样率必须至少是信号带宽的 2 倍。考虑一个由多个信号（可能是简单的正弦波）组成的信号，其中最高频率的正弦波为 50 Hz。为了捕捉所有相关信息，采样率必须在 100 Hz 或更高。有关更深入的解释，请参阅下文。压缩感知的目的是允许我们在仍然能够完美（或近乎完美）地恢复信号的情况下，采集比之前所需样本量更少的样本。

http://en.wikipedia.org/wiki/Compressed_sensing

http://en.wikipedia.org/wiki/Nyquist_rate

背景

当我第一次听说“压缩感知”时，我觉得它非常有趣，但绝对没有一个扎实的理解。我写这篇文章不仅是为了启发他人了解 CS，也是为了加深我自己对它的理解。以下是一些可能需要预先了解的领域：

线性代数/矩阵代数

由于我们将使用 Matlab，线性代数总是很有用的，特别是矩阵的处理方式。了解 ' 与 . 以及 * 与 .* 之间的区别很重要（查阅 Matlab 本身的 复习页面）。

信号处理基础

同样有帮助的是一些信号处理知识，例如傅里叶变换和奈奎斯特定理（见导言）的基本概念。我们将要观察的数据可能代表简单的时域信号，但理解什么是变换很重要，因为这是 CS 的核心。

压缩技术 (JPEG)

这稍微有点高级，可能不是一开始就必须理解的。但是，如果你曾经研究过计算机上的图像和文件是如何变得如此小的，你可能已经遇到过压缩。对于图像压缩，JPEG 是一个旧的常用标准。JPEG 使用一种称为离散余弦变换 (DCT) 的技术。后来创建了一个不同的版本，称为 JPEG2000，它使用小波变换。这些变换是我们稍后将要介绍的一些更高级的方法。

稀疏性

稀疏性是压缩感知中的一个常用术语。在数学意义上，一个**稀疏**的数据集合具有少量的非零值。简单来说，大多数值是 0，只有少数包含有意义的数据。请查看关于稀疏矩阵的条目。

压缩感知背后的基本数学原理

现在，让我们深入探讨 CS 背后的数学原理。我将尽量使用我在文献中看到的一些常用符号，这样如果你决定进一步阅读，就不会那么吃力了。

想象我们有一个信号 *x*，它代表一个长度为 N 的时域信号。在理想情况下，我们将简单地采集 N 个测量值，然后就完成了！然而，今天我们只能进行 K 次测量（K < N）。

我们使用一个向量 *y* 来表示我们的 K 次测量。需要注意的是，**一次测量并不一定对应于一个输入值**。起初这可能令人困惑，但现实生活中的测量并非总是以单一、线性的方式进行。

为了得到 *y*，我们使用 *A*。简单地说：*y* = *Ax*。

*A* 是**传感矩阵**，它允许我们通过（随机测量¹、变换或两者的组合）从 *x* 得到 *y*。

总结一下，我们有以下内容：

• *x* - N×1 向量，表示原始信号（必须是某种可压缩的）
• *y* - K×1 向量，表示输出值
• *A* - K×N 传感矩阵（如何从输入到输出）

 

练习 1

假设 *x* 的长度为 N=8，并包含随机值。我们只想保留其中的 3 个（假设是第 1、2 和第 5 个）。在 *y* = *Ax* 的情况下，线性代数看起来会是这样。

在图 1 中，我们显然保留了原始值的 3 个样本，并丢弃了其余 5 个。如果 *x* 是一个稀疏信号（只有少数有用的值），那么我们就不能随机地保留几个值，否则我们可能会丢失我们需要的那些！这是因为我们事先不知道有用值的确切位置。

_{图 2：如果我们从稀疏信号中随机选择值，我们就会丢失一些信息！}

我们需要一种不同的方法来避免丢失我们的非零值！幸运的是，如果我们使用高斯随机变量作为我们的传感矩阵 *A*，我们就不会丢失数据。目前它看起来不一样，但它仍然编码在我们的输出中。

_{图 3：以不同的方式获得仅 3 个输出值 - 随机测量。}

 

在图 3 所示的系统中，我们有一些已知量、一些未知量和一些假设：

已知

• 3×8 传感矩阵 *A*（它可能是随机生成的，但*我们*生成了它）
• 3×1 的输出向量 *y*

未知

• 原始值 *x*

假设

• 8×1 向量 *x* 是稀疏的（如果它不稀疏，我们所做的就无效）

有了以上知识，我们如何确定 *x* 是什么？这就是压缩感知恢复方法发挥作用的地方。一种流行的方法是 ℓ₁-MAGIC 算法套件 [1]，它以 L¹ 范数最小化² 为中心，超出了本文的范围。这是一个迭代过程，依赖于多个参数，但上述条件是满足该算法要求的良好开端。

Matlab 代码

请注意，您必须拥有 [1] 中的 l1magic 代码文件夹才能运行这些示例。我建议下载代码并将其放在您用于以下示例的同一个目录中。（遗憾的是，这里没有 Matlab 的代码模板，所以您只能看到单色文本）。

重要提示

您需要更改文件 "l1eq_pd.m" 的第 141 行以避免错误。这一行最初是这样的：

[dv,hcond] = linsolve(H11p, w1p, opts);

我们希望将其更改为以下内容：

[dv,hcond] = linsolve(H11p, w1p);

我建议您只需注释掉第 141 行，并在其下方添加更正后的行。

示例 1：稀疏时域信号

% Initialize constants and variables
rng(0);                 % set RNG seed
N = 256;                % length of signal
P = 5;                  % number of non-zero peaks
K = 64;                 % number of measurements to take (N < L)
x = zeros(N,1);         % original signal (P-sparse)

% Generate signal with P randomly spread values
peaks = randperm(N);
peaks = peaks(1:P);
x(peaks) = randn(1, P);
amp = 1.2*max(abs(x));
figure; subplot(3,1,1); plot(x); title('Original signal'); xlim([1 N]); ylim([-amp amp]);

% Obtain K measurements
A = randn(K, N);
y = A*x;
subplot(3,1,2); plot(y); title('K measured values'); xlim([1 K]);

% Perform Compressed Sensing recovery
x0 = A.'*y;
xp = l1eq_pd(x0, A, [], y);
subplot(3,1,3); plot(real(xp)); title('Recovered signal'); xlim([1 N]); ylim([-amp amp]);

原始信号有 256 个值，但我们完美地恢复了它，即使我们只采集了 64 个随机测量值！

正如您所见，压缩感知在这个简单的例子中是有效的。我建议您尝试一下，看看结果是否会改变。例如，如果我们采集的测量值更少会怎么样？尝试将 K 减小到 32 或 16，看看会发生什么。

更多压缩感知数学原理

有时我们感兴趣的信号在其测量的域中并不是稀疏的。例如，考虑一个由少量正弦波组成的信号。在时域中，函数 sin(2πt) 不是稀疏的。但是，如果我们观察它在频域中的情况，它完全由两个峰表示³。

为了在稀疏域中执行恢复，我们必须明确地这样实现。为了做到这一点，理解对于 CS 至关重要的两个标准正交基（或正交基）非常重要。

理解 CS 中的正交基

Ψ - **表示基**

如前所述，为了执行压缩感知，我们的信号必须是稀疏的。表示基 Ψ 是将信号转换到稀疏域的数学方法。

在示例 1 中，时域信号已经在其稀疏域中。因此，Ψ 只是单位矩阵 **I**_N，可以从代数步骤中排除。然而，在下面的示例 2 中，我们将有一个在不同域（在这种情况下是频域）中稀疏的时域信号。因此，Ψ 将是那个将时域信号转换为频域的基，即DFT 矩阵。

如 [2] 所述，**Ψ** 有许多可能性，选择正确的取决于具体应用。对于频域稀疏信号，DFT 是适用的。然而，对于 2D 图像，许多小波变换基之一可能适用。

Φ - **传感基**

传感基 Φ 表示我们从中提取信号值的域。在示例 1 中，这只是具有随机高斯条目的矩阵。在下面的示例 2 中，它将是脉冲。正如您所看到的，这描述了用于获取较小（压缩）数据量的方法。在示例 1 和 2 中，此操作分两个步骤进行：

1. 创建一个具有 N 个值的随机排列。（示例 2 中有两个 `randperm`，请务必找到正确的那个）。
2. 获取前 K 个值，确保我们随机获得了 K 个 Φ 的子集。

附加信息

为了简化起见，我们一直将结果 K×N 矩阵（例如，示例 1 中的 *A*）视为传感基。然而，为了精确起见，传感基实际上是完整的 N×N 矩阵，并且使用另一个矩阵 *R* 从 Φ 中提取 K 个条目。提取矩阵是通过简单地保留 I_N 中的 K 行获得的。

因此，传感矩阵的完整形式为 *A* = *R*ΦΨ。

这对我们的代码造成的影响是，恢复的数据将不在我们原来的域中。它在压缩域中，因此必须通过乘以 Ψ^-1 来“解压缩”。这将数据返回到所需的域。

Matlab 代码（续）

示例 2：频域稀疏信号

代码前的几点说明

您可以通过在频域中生成一组固定数量的脉冲而不是先生成时域信号来轻松执行相同的示例。下面方法的好处在于，它使得使用实值时域信号更容易，从而更准确地表示现实生活中的信号。

因为在这种情况下 Ψ 是 DFT 矩阵，它的逆很容易获得，即 conj(Ψ)/N （请参阅 dftmtx 的帮助页面）。

如果您没有 Matlab 的信号处理工具箱，您可以将 `dftmtx` 替换为类似的算法（例如这个），甚至为了练习自己编写算法！（类似于这里，但没有 *N*^-1/2 缩放因子）

% Initialize constants and variables
rng(0);                 % set RNG seed
N = 256;                % length of signal
P = 2;                  % number of sinusoids
K = 64;                 % number of measurements to take (N < L)

% Generate signal with P randomly spread sinosoids
% Note that a real-valued sinusoid has two peaks in the frequency domain
freq = randperm(N/2)-1;
freq = freq(1:P).';
n = 0:N-1;
x = sum(sin(2*pi*freq/N*n).', 2);

% Orthonormal basis matrix
Psi = dftmtx(N);
Psi_inv = conj(Psi)/N;
X = Psi*x;              % FFT of x(t)

% Plot signals
amp = 1.2*max(abs(x));
figure; subplot(5,1,1); plot(x); xlim([1 N]); ylim([-amp amp]);
title('$\mathbf{x(t)}$', 'Interpreter', 'latex')
subplot(5,1,2); plot(abs(X)); xlim([1 N]);
title('$|\mathbf{X(f)}|$', 'Interpreter', 'latex');

% Obtain K measurements
x_m = zeros(N,1);
q = randperm(N);
q = q(1:K);
x_m(q) = x(q);
subplot(5,1,3); plot(real(x_m)); xlim([1 N]);
title('Measured samples of $\mathbf{x(t)}$', 'Interpreter', 'latex');

A = Psi_inv(q, :);      % sensing matrix
y = A*X;                % measured values

% Perform Compressed Sensing recovery
x0 = A.'*y;
X_hat = l1eq_pd(x0, A, [], y, 1e-5);

subplot(5,1,4); plot(abs(X_hat)); xlim([1 N]);
title('$|\mathbf{\hat{X}(f)}|$', 'Interpreter', 'latex');

x_hat = real(Psi_inv*X_hat);    % IFFT of X_hat(f)

subplot(5,1,5); plot(x_hat); xlim([1 N]);  ylim([-amp amp]);
title('$\mathbf{\hat{x}(t)}$', 'Interpreter', 'latex');

同样，原始时域信号有 256 个值。但是，通过只从中提取 64 个单独的值（25%），我们能够恢复所有 256 个原始值！

这与示例 1 不同之处在于，感兴趣的信号不是在时间域中稀疏，而是在*频域*中稀疏。之所以有效，是因为传感基（脉冲）和表示基（DFT）彼此正交。一如既往，我鼓励您通过更改参数（如 N、K 和 P）来进行实验。

关注点

此处包含的示例已使用预设值对随机数生成器 (RNG) 进行了播种。这是为了确保您获得与所示相同的​​结果。您可以随意更改种子值（当前为 0）或完全注释掉该行，以便每次使用不同的随机变量进行测试。

Matlab 帮助 (F1) 是您的好帮手，因此当您需要帮助理解某个函数时，请务必查阅它。

别忘了仔细检查您的转置运算符。处理实值数据时，' 和 . 是相同的。
[1 2 3] 的共轭转置是 [1 2 3]^T，这与转置*相同*。
然而，[1i 2i 3i] 的共轭转置变为 [-1i -2i -3i]^T，这与转置*不相同*。

脚注

¹ “随机测量”实际上不一定需要是随机的。它们可以是预定的序列，如伪随机码、噪声小波、二进制码等。然而，目前这超出了本文的范围。

² https://mathworld.net.cn/L1-Norm.html

³ 复指数正弦波在频域中是一个峰，但实值指数正弦波等同于频率为 f 和 -f 的两个复指数正弦波，因此在频域中有两个峰。

这也取决于采样频率。如果正弦波的频率不是 1/N 的整数倍，那么正弦波在频域中会展宽，而不是一个尖锐的峰（请参阅相关的Matlab 文章）。

历史

V1：导言、背景、数学、一维时域稀疏示例、一维频域稀疏示例。

未来可能的添加内容

• 二维示例。
• 带噪声的恢复

参考文献

1. E. Candès 和 J. Romberg，ℓ₁-MAGIC
   网址：http://statweb.stanford.edu/~candes/l1magic/
2. Candes, E.J.; Wakin, M.B., "An Introduction To Compressive Sampling," Signal Processing Magazine, IEEE, vol.25, no.2, pp.21,30, March 2008. doi: 10.1109/MSP.2007.914731
   网址： http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=4472240&isnumber=4472102
3. Matlab 帮助文件。Mathworks, Cambridge MA. 访问于 2014 年。
   网址：http://www.mathworks.com/help/matlab/