CN103250352A - 用于测量和恢复稀疏信号的方法和装置 - Google Patents

Abstract

本发明提供两种新压缩感测技术。在第一技术解决方案中，提出基于排列的多维感测矩阵和具有最大似然(ML)局部检测的迭代恢复算法，其能够完全利用稀疏信号的数字性质。在第二技术解决方案中，提出包含基于排列的多维测量矩阵的稀疏测量矩阵和完全利用测量符号的特征以设计每个迭代中的简单局部恢复的迭代恢复算法。第二技术解决方案可以同时实现线性解码复杂度并且经验地实现概略长度下界。

Description

用于测量和恢复稀疏信号的方法和装置

技术领域

本申请涉及压缩感测技术。

背景技术

压缩感测(CS)是一种近来开发的技术。考虑自然和人为信号的大部分具有稀疏或者接近稀疏性质这样的事实，压缩感测技术可以发现在许多不同领域、比如压缩成像、压缩采样、信号处理、数据流计算和组合组测试中的应用。压缩感测的基本思想在于可以从长度为M的线性测量y＝Ax准确恢复长度为N的稀疏信号x(信号如果它包含比非零元素多得多的零元素则称为稀疏)，其中A是M×N测量矩阵，M＜＜N。

可以通过使说明测量矢量的‖x‖₀最小化来执行重建。由于这一最小化问题是NP难度，所以已经考察次优算法。计算上可行的稀疏信号恢复算法的主要类别包括通过经常选择p为1的l_p最小化问题近似最小化l₀问题并且使用这一凸优化来求解这一问题的凸松弛；通过相继地标识产生最大质量提高的一个或者多个分量来迭代地精化稀疏解的匹配追逐；以及假设有利于信号矢量稀疏性的先验分布并且使用最大后验估计器以结合观测的贝叶斯框架。尽管它们在实践中有相对良好性能，但是它们最适合于具有连续值的信号。对于具有数字值的稀疏信号，例如在处理单色图像时，这些算法有欠充分，因为它们不能利用源的数字性质，该数字性质如果被恰当利用则可以大量提高恢复准确度。

因此需要一种可以完全利用信号的数字性质的新压缩感测技术。

此外，在几乎所有应用中，优选的是测量矩阵A是稀疏的、即它在每列中包含比非零条目多得多的零条目。稀疏测量矩阵的优点包括在编码和解码二者中的低计算复杂度、对信号的增量更新容易和低存储要求等。许多研究已经致力于利用稀疏测量矩阵的CS，但是它们中的多数无法同时实现线性解码复杂度和性能限界。现有算法的典型例子包括匹配追逐和凸优化。匹配追逐型算法可以用线性恢复复杂度渐近地实现概略长度的下界。然而数值结构已经示出在这一类算法中需要的经验概略长度总是比渐进限界高得多。在另一方面，凸优化型算法可以渐进地和经验地实现概略长度的下界，这在实践中表现在测量数目方面的优点。例如用非零元素数目K＝50和信号长度N＝20000，示出匹配追逐需要约2000个测量而凸优化仅需约450个测量。凸优化型算法的一大缺点是它们的随着信号长度N以多项式阶O(N³)增长的更高恢复复杂度。

因此需要一种可以用稀疏测量矩阵同时实现线性解码复杂度并且经验地实现概略长度下界的新压缩感测技术。

发明内容

为了更好地解决上述两个问题，提供两个技术解决方案，一个是提供一种可以完全利用稀疏信号的数字性质的新压缩感测技术，而另一个是提供一种可以用稀疏测量矩阵同时实现线性解码复杂度并且经验地实现概略长度下界(lower bound of sketch length)的新压缩感测技术。

基于这一点，在本发明的第一方面中，提供一种用于处理数字稀疏信号的方法。该方法包括以下步骤：用M×N测量矩阵A对长度为N的稀疏信号矢量x＝{x_i}执行线性测量以获得长度为M的测量矢量y，其中测量矩阵A由下式表示：

$A=\left[\begin{array}{c}{A}_1{\mathrm{\Π}}_1\\ \·\\ \·\\ \·\\ A_d{\mathrm{\Π}}_d\\ \·\\ \·\\ \·\\ A_D{\mathrm{\Π}}_D\end{array}\right],$

其中A_d∏_d是测量矩阵A的子矩阵，其中d＝1～D，∏_d是N×N随机排列矩阵，并且A_d是下式表示的J×N矩阵：

其中M＜＜N，J*L＝N，并且稀疏信号矢量x中的每个条目x_i被取自于有限集合Q＝{X₀＝0，X₁，...，X_Q-1}而X_q为非零数，其中q＝1～Q-1，并且Q是集合的大小，其中根据相同排列矩阵生成的测量符号称为一个维度，并且D是总维度数目。

有利地，该方法还可以包括以下步骤：使用D个最大似然检测器以针对D个维度分别执行最大似然检测，其中第d个最大似然检测器用于针对第d个维度执行最大似然检测；将上述步骤重复多个迭代直至满足预定条件；并且基于第D个最大似然检测器在最后迭代中的输出来估计稀疏信号矢量中的源符号。

在本发明的第一方面中的测量矩阵A允许在每个维度中的简单最大似然检测，该最大似然检测完全利用稀疏信号的数字性质并且针对每个维度提供计算上可行的局部优化检测。测量矩阵的多维结构实现在维度之间的迭代信息交换以获得接近全局最优的估计结果。

另外，在本发明的第一方面中的测量矩阵A是稀疏的、即它在每列中包含比非零条目多得多的零条目。测量矩阵的稀疏具有若干有吸引力的性质、比如在编码和恢复二者中的低计算复杂度、对信号的增量更新容易和低存储要求等。这些优点使第一方面中的技术解决方案成为用于利用稀疏数字信号的压缩感测的潜在和实用解决方案。

在本发明的第二方面中，提供一种用于处理模拟稀疏信号的方法。该方法包括以下步骤：用M×N测量矩阵A对长度为N的K稀疏信号矢量x执行线性测量以获得长度为M的测量矢量y，其中测量矩阵A由下式表示：

$A=\left[\begin{array}{c}{A}_1{\mathrm{\Π}}_1\\ \·\\ \·\\ \·\\ A_d{\mathrm{\Π}}_d\\ \·\\ \·\\ \·\\ A_D{\mathrm{\Π}}_D\end{array}\right],$

其中A_d∏_d是测量矩阵A的子矩阵，其中d＝1～D，∏_d是N×N随机排列矩阵，并且A_d是下式表示的J×N矩阵：

其中K＜＜N，M＜＜N，J*L＝N，并且稀疏信号矢量

其中根据相同排列矩阵生成的测量符号称为一个维度，并且D是总维度数目。

有利地，该方法还可以包括以下步骤：

i.针对

中的每个元素

其中d＝1～D，j＝1～J，判断是否等于0，其中初始化

为

并且

是

的第d个维度中的第j个元素；

ii.如果

等于0，则设置

其中l＝1～L；其中∏_d(i)是x_i在第d个排列版本中的索引，并且

是∏_d(i)的逆运算；

并且该方法还包括以下步骤：

u.针对稀疏信号矢量x中的每个源符号x_i，其中i＝1～N，判断是否满足以下等式：

$\frac{{\tilde{y}}_{j_i^{\left(d_m\right)}}^{\left(d_m\right)}}{a_{j_i^{\left(d_m\right)},l_i^{\left(d_m\right)}}^{\left(d_m\right)}}=\frac{{\tilde{y}}_{j_i^{\left(d_n\right)}}^{\left(d_n\right)}}{a_{j_i^{\left(d_n\right)},l_i^{\left(d_n\right)}}^{\left(d_n\right)}},\∀d_m=1~D,d_n=1~D,d_m\≠d_n;$

v.如果针对一对(d_m，d_n)满足等式，则设置

其中 $j_i^{\left(d_m\right)}=\left|({\mathrm{\Π}}_{d_m}\left(i\right)-1/L)\right|$ 并且 $l_i^{\left(d_m\right)}=\mathrm{mod}({\mathrm{\Π}}_{d_m}\left(i\right)-1,L)+1;$

并且在上述步骤之后，该方法还包括以下步骤：

p.按照以下公式通过从测量矢量y减去恢复的稀疏信号来更新

$\tilde{y}=y-A\tilde{x},$

其中

代表其中未恢复的符号设置为0的恢复的稀疏信号矢量；

q.将上述步骤重复多个迭代直至满足预定条件。

在本发明的第二方面中的测量矩阵A的特殊结构造成测量符号的如后文将说明的一些感兴趣的特征，可以利用这些特征以设计每个迭代中的很简单的恢复算法。迭代过程用来从更易于恢复的信号开始并且取消已经恢复的符号的贡献以有助于恢复其它符号来逐步恢复源信号。通过重复这样的检测和取消操作，可以获得接近全局最优解。需要的复杂度仅随着源信号长度N线性地增长。多维结构和随机排列矩阵确保每个测量向所有源符号在统计上(直接或者间接)提供对于提出的技术的良好性能而言关键的一些有用信息。

另外，在本发明的第二方面中的技术解决方案可以同时实现经验概略长度下界和线性复杂度。良好经验性能和低复杂度使它成为用于利用稀疏测量矩阵的压缩感测的良好和实用备选解决方案。

此外，根据本发明的一个实施例，提供一种用于测量和恢复稀疏信号的装置，该装置包括：测量装置，用于用M×N测量矩阵A对长度为N的稀疏信号矢量x＝{x_i}执行线性测量以获得长度为M的测量矢量y，其中测量矩阵A由下式表示：

$A=\left[\begin{array}{c}{A}_1{\mathrm{\Π}}_1\\ \·\\ \·\\ \·\\ A_d{\mathrm{\Π}}_d\\ \·\\ \·\\ \·\\ A_D{\mathrm{\Π}}_D\end{array}\right],$

其中A_d∏_d是测量矩阵A的子矩阵，其中d＝1～D，∏_d是N×N随机排列矩阵，并且A_d是下式表示的J×N矩阵：

其中M＜＜N，J*L＝N，根据相同排列矩阵生成的测量符号称为一个维度，并且D是总维度数目；以及

恢复装置，用于从长度为M的测量矢量y恢复长度为N的稀疏信号矢量x＝{x_i}。

附图说明

参照附图更具体地并且通过例子说明本发明，在附图中：

图1示出根据一个实施例的处理数字稀疏信号的流程图；

图2示出图1中使用的迭代算法的原理；

图3a示出在p₁＝0.1时在图1的实施例与两种常规技术之间的速率-失真性能比较；

图3b示出在p₁＝0.05时在图1的实施例与两种常规技术之间的速率-失真性能比较；

图4示出根据另一实施例的处理模拟稀疏信号的流程图；

图5示出在图4的实施例的测量数目M与Klog₂(N/K)之间的关系；

图6a示出在p₁＝0.1时在图4的实施例的恢复失败概率与α之间的关系；并且

图6b示出在p₁＝0.05时在图4的实施例的恢复失败概率与α之间的关系。

贯穿上述附图，将理解相似标号指代相似、类似或者对应特征或者功能。

具体实施方式

下文将分别更具体描述本发明的上述第一方面的技术解决方案和本发明的上述第二方面的技术解决方案。

图1示出根据一个实施例的处理数字稀疏信号的流程图。

在步骤S11中，用M×N测量矩阵A对长度为N的数字稀疏信号矢量x执行线性测量(即编码)以获得长度为M的测量矢量y，其中M＜＜N。

考虑具有独立和相同分布(i.i.d.)的条目的数字稀疏信号矢量x＝{x_i}。稀疏信号矢量x中的每个条目x_i的值被取自于有限集合Q＝{X₀＝0，X₁，...，X_Q-1}而X_q为非零数，其中q＝1～Q-1，并且Q是集合的大小。假设每个条目x_i具有为“0”的概率p₀和为X_q的概率p_q(q＝1～Q-1)。由于x是稀疏的，所以

基于压缩感测理论，可以从测量矢量y恢复x如下：

y＝Ax+n    (1)

其中n是具有零均值和E(|n|²)≤Mσ²的长度为M的噪声矢量。

在这一实施例中，设计测量矩阵A如下：

$A=\left[\begin{array}{c}{A}_1{\mathrm{\Π}}_1\\ \·\\ \·\\ \·\\ A_d{\mathrm{\Π}}_d\\ \·\\ \·\\ \·\\ A_D{\mathrm{\Π}}_D\end{array}\right]---\left(2\right)$

其中A_d∏_d(d＝1～D)是测量矩阵A的子矩阵，∏_d是N×N随机排列矩阵，并且A_d是设计如下的J×N矩阵：

其中J*L＝N，根据相同排列矩阵生成的测量符号称为一个维度，并且D是总维度数目。针对D个维度独立和随机生成排列矩阵{∏_d}。{A_d}中的非零条目是高斯随机变量。

可以实施编码过程如下。独立排列数字稀疏信号矢量x中的源符号D次。将每个排列版本划分成各自包含L个符号的J＝N/L组。在第d个排列版本的第j组中的符号被A_d的第j行中的对应非零条目加权、然后被线性地叠加以生成这一维度中的第j个测量符号。测量总数M由组长度L、维度数目D和信号长度N确定为M＝ND/L。

在编码之后，应用量化过程以数字化测量符号。量化误差由公式(1)中的噪声矢量n代表。令S表示量化级数并且p_quan(s)是测量符号被量化成S级的概率。为了代表一个量化测量而需要的位数然后是：

$b=\underset{s=1}{\overset{S}{\mathrm{\Σ}}}{p}_{\mathrm{quan}}\left(S\right)\·{\log }_2(1/p_{\mathrm{quan}}\left(s\right))---\left(4\right)$

因此，为了代表x而需要的总位数是B＝bM＝bND/L，并且每个源符号需要的平均位数是：

$\mathrm{\η}=\frac{B}{N}=\mathrm{bD}/L---\left(5\right)$

可以经由选择D、L和S来调整每个源符号的平均位数。

公式(3)中的{A_d}的结构实现将最大似然检测(ML)用于每个维度中的每组L个符号。通过为L选择小值，可以控制ML检测的计算复杂度。ML检测可以完全利用源信号的数字性质并且提供用于每个维度的局部最优解。多维结构允许在维度之间的迭代信息交换以实现接近最优的估计。在不同维度中的独立随机排列矩阵确保每个测量向数字稀疏信号矢量x中的所有符号在统计上(直接或者间接)提供一些有用信息，因为对一个维度中的不同测量有贡献的符号可以被分组在一起并且对其它维度中的相同测量有贡献。这意味着一个测量向它的关联L个符号提供的信息可以帮助检测其它符号，如果这些L个符号与其它维度中的其它符号分组。在迭代恢复算法中完全利用这样的性质以大量提高提出的技术的抗噪声能力。

基于这一点，在步骤S12中，D个最大似然检测器用来针对D个维度分别执行最大似然检测。也就是说，如图2中所示，每个最大似然检测器负责检测一个维度。图2示出迭代算法的原理，其中“DET-d”是第d个维度的局部检测器，“T”表示一个迭代的延迟，并且“/”表示除法运算。

在每个局部检测器内，基于每组L个符号的噪声测量符号和先验信息对每组L个符号执行ML检测。每个检测器的输出是数字稀疏信号矢量x中的源符号的后验信息，该后验信息用来在下一迭代中精化其它维度中的局部检测。

定义图2中涉及到的变量如下。

在第d个维度中x_i为X_q(q＝0～Q-1)的先验概率。

p^(d)(x_i＝X_q)：在第d个维度中x_i为X_q(q＝0～Q-1)的后验概率。

e^(d)(x_i＝X_q)：在第d个维度中x_i为X_q(q＝0～Q-1)的外赋概率。

在第一迭代中初始化先验概率为

其中i＝1～N并且q＝0～Q-1，并且在第一迭代中初始化外赋概率为e^(d)(x_i＝X_q)＝1，其中

考虑第d个维度中的由

表示的特定测量j。令

是对

有贡献的L个源符号的索引。

的值由排列矩阵∏_d确定。使用ML检测来计算这些L个符号的后验概率为：

$p^{\left(d\right)}(x_{i_{j,l}^{\left(d\right)}}=X_q)$

$=\frac{1}{\sqrt{{2\mathrm{\π\σ}}^2}}\underset{c_l=X_q}{\underset{c\∈Q^L}{\mathrm{\Σ}}}\exp \{-\frac{{|y_j^{\left(d\right)}-\underset{l^{\′}=1}{\overset{L}{\mathrm{\Σ}}}{a}_{j,l^{\′}}^{\left(d\right)}{c}_{l^{\′}}|}^2\·\underset{l^{\′}=1}{\overset{L}{\mathrm{\Π}}}{\tilde{p}}^{\left(d\right)}(x_{i_{j,l^{\′}}^{\left(d\right)}}=c_{l^{\′}})}{{2\mathrm{\σ}}^2}\},$ 其中l＝1～L  (6)

公式(6)中的求和是针对所有可能矢量c∈Q^L而第l个元素被固定成X_q。通过从后验概率提取先验概率来计算外赋概率，

$e^{\left(d\right)}(x_i=X_q)=p^{\left(d\right)}(x_i=X_q)/{\tilde{p}}^{\left(d\right)}(x_i=X_q)---\left(7\right)$

在第d个维度中生成的后验概率用来更新第(mod(d，D)+1)个维度中的先验概率为：

${\tilde{p}}^{(\mathrm{mod}(d,D)+1)}(x_i=X_q)=p^{\left(d\right)}(x_i=X_q)/e^{(\mathrm{mod}(d,D)+1)}(x_i=X_q)---\left(8\right)$

注意在先前迭代中在第(mod(d，D)+1)个维度中生成e^{(mod(d，D)+1)}(x_i＝X_q)，因此可以防止它再次循环回到第(mod(d，D)+1)个维度，这是迭代检测的基本规则。这由公式(8)中的除法运算实现。

然后在步骤S13中，将上述过程重复多个迭代直至满足预定条件。

预定条件可以例如包括以下任一项：

-已经达到固定迭代数目；

-在两个相继迭代中在第D个最大似然检测器中生成的后验概率之差在预定义阈值以下。

最后在步骤S14中，基于第D个最大似然检测器在最后迭代中的输出来估计数字稀疏信号矢量x中的源符号{x_i}。

例如针对源符号{x_i}进行硬判决为：

${\hat{x}}_i=X_{\widehat{\mathrm{qi}}}---\left(9\right)$

其中

${\hat{q}}_i=\underset{q}{\arg \max }\left(\right\{p^{\left(D\right)}(x_i=X_q)|q=0~Q-1\})---\left(10\right)$

上文提出的迭代算法的复杂度在这一实施例中是O(NDQ^L/L)，这比最优解O(Q^N)的复杂度低得多，因为L比N小得多并且未随N增长。注意更小L将针对固定D造成更大测量数目M。然后通过调整D和L的值，可以实现在计算复杂度与测量数目之间的不同权衡。

在下文中，使用数值结果来举例说明上文提出的技术解决方案的性能。考虑数字稀疏信号矢量x中的二进制源信号而x_i∈{0，1}。x的条目是i.i.d.变量而为“0”(“1”)的概率p₀(p₁)且p₀＞＞p₁。根据公式(2)和(3)生成测量矩阵A而非零条目为具有零均值和单位方差的i.i.d.高斯随机变量。归一化测量矩阵A，从而每行具有单位范数。用以下量化规则将量化级固定成S＝5，其中Quan(x)代表x在量化之后的值。

注意就数字稀疏信号矢量x而言，使用在公式(2)和(3)中定义的测量矩阵A来生成的线性测量符号(在量化之前)具有为“0”的大概率。如果零测量符号出现，则可以用高准确度推断与它关联的所有L个源符号为零，这将很有助于在迭代恢复算法中检测其它源符号。因此，用于公式(20)中的“0”的具体量化级用来区分它与其它值。使用MonteCarlo仿真以获得概率{p_quan(s)，s＝1～5}并且根据公式(5)计算每个源符号需要的平均位数η。

图3a和3b分别示出在上文提出的技术解决方案与基于随机高斯感测矩阵、Lasso算法和贝叶斯框架的两种常规技术之间的速率-失真性能比较。也示出用于无损恢复的每个源符号需要的最大位数用于参考。

这里设置p₁＝0.1和0.05并且将迭代数目固定于5。在表1中列举概率{p_quan(s)}。

表1.在量化之后的概率{p_quan(s)}。

从表1可以计算每个源符号需要的位数η。将D从1调整到4并且将L从10调整到12以获得从0.14到1的不同η值。测量ξ表示的失真如下，该失真等于由N归一化的在中的不正确条目平均数目。

$\mathrm{\ξ}=\frac{1}{N}E\left({\left|\right|\hat{x}-x\left|\right|}_2^2\right)---\left(12\right)$

作为参考，在图3中针对经由熵编码的无损压缩示出每个源符号需要的最小位数。也包括两种常规方式的速率-失真性能，其基于随机高斯感测矩阵并且分别运用凸松弛(细节请参见[R.Gribonval和M.Nielsen，“Highly sparse representations from dictionaries are unique andindependent of the sparseness measure，”Aalborg Univ.，Aalborg，Denmark，Tech.Rep.，Oct.2003.]和[J.A.Tropp，“Just relax：Convex programmingmethods for identifying sparse signals in noise，”IEEE Trans.Inf.Theory，vol.52，no.3，pp.1030-1051，Mar.2006.])以及贝叶斯框架(细节请参见[S.Ji，Y.Xue and L.Carin，“Bayesian compressive sensing，”IEEE Trans.Signal Processing，vol.56，no.6，pp.2346-2356，June 2008.]和[M.E.Tipping，“Sparse Bayesian learning and the relevance vector machine，”Journal of Machine Learning Research，vol.1，pp.211-244，Sept.2001.])用于信号恢复。对于凸松弛，使用Lasso算法(细节请参见M.J.Wainwright，“Sharp thresholds for high-dimensional and noisy sparsityrecovery using l1-constrained quadratic programming(Lasso)，”IEEE Trans.on Inform.Theory，vol.55，pp.2183-2202，May 2009.]和[R.Tibshirani，“Regression shrinkage and selection via the lasso，”Journal of the RoyalStatistical Society，Series B，pp.267-288，1996.])。对于基于贝叶斯框架的算法，使用论文{S.Ji，Y.Xue and L.Carin，“Bayesian compressivesensing，”IEEE Trans.Signal Processing，vol.56，no.6，pp.2346-2356，June2008.’}中介绍的技术。两种常规方式除了去除专用于“0”的量化级之外使用与公式(11)中相似的量化规则，因为在这一情况下，每个测量很不可能为“0”。在表1中也列举用于常规方式的概率{p_quan(s)，s＝1～5}。从图3观测到提出的技术可以实现比基于随机感测矩阵和凸松弛/贝叶斯框架的常规方式好得多的速率-失真性能。

图4示出根据另一实施例的处理模拟稀疏信号的流程图。

在步骤S41中，用M×N测量矩阵A对长度为N的K稀疏信号矢量x执行线性测量(即编码)以获得长度为M的测量矢量y，其中K＜＜N，M＜＜N，并且K稀疏信号矢量

测量矢量

根据压缩感测理论，可以从测量矢量y恢复x如下：

y＝Ax    (13)

在这一实施例中，设计测量矩阵A如下：

$A=\left[\begin{array}{c}{A}_1{\mathrm{\Π}}_1\\ \·\\ \·\\ \·\\ A_d{\mathrm{\Π}}_d\\ \·\\ \·\\ \·\\ A_D{\mathrm{\Π}}_D\end{array}\right]---\left(14\right)$

其中A_d∏_d(d＝1～D)是测量矩阵A的子矩阵，∏_d是N×N随机排列矩阵，并且A_d是设计如下的J×N矩阵：

其中J*L＝N，根据相同排列矩阵生成的测量符号称为一个维度，并且D是总维度数目。显然，以这一方式生成的测量矩阵A是稀疏的而每列仅有D个非零条目。

可以实施编码操作如下。独立排列数字稀疏信号矢量x中的源符号D次。将每个排列版本划分成各自包含L个符号的J＝N/L组。在第d个排列版本的第j组中的符号被A_d的第j行中的对应非零条目加权、然后被线性地叠加以生成这一维度中的第j个测量符号。

表示第d个维度中的第j个测量符号，该测量符号具有：

$y_j^{\left(d\right)}=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{a}_{j,l}^{\left(d\right)}{\hat{x}}_{(j-1)L+1}^{\left(d\right)}---\left(16\right)$

其中

是x的第d个排列版本。令∏_d(i)是x_i在第d个排列版本中的索引，并且

是逆运算，它们具有：；

${\hat{x}}_{{\mathrm{\Π}}_d\left(i\right)}^{\left(d\right)}=x_i$ 和 ${\hat{x}}_i^{\left(d\right)}=x_{{\mathrm{\Π}}_d^{-1}\left(i\right)}---\left(17\right)$

稀疏信号矢量x中的每个源符号x_i与各自在一个维度中的D个测量关联。总测量数目(即概略长度)M由组长度L、维度数目D和源信号长度N确定为M＝ND/L。

从公式(16)可见测量矩阵A在每列中仅有D个非零元素。也注意x的仅非零元素造成编码中的加法和乘法运算。因此，提出的技术的编码复杂度约为DK个乘法和DK个加法。

公式(15)中的{A_d}的块对角结构和D个维度中的随机排列矩阵造成测量符号的以下特征，这在信号恢复中很有用。为了便于讨论，定义测量符号的度为与它关联的非零源符号的数目。从(16)可见每个测量符号

与仅L个源符号

关联，因此易于验证以下两个问题具有相对大的发生概率。

问题-1：测量符号具有0度。

问题-2：对于非零源符号，它的关联测量中的至少两个测量具有1度(假设D≥2)。

可以计算上述两个问题的概率如下。令p₁＝K/N并且p₀＝1-p₁。p₁和p₀分别代表源符号为非零和零的概率。然后可以近似地计算问题1和问题2的概率为：

$P_{\mathrm{issue}-1}\≈p_0^L={\left(\frac{N-K}{N}\right)}^L---\left(18\right)$

并且

$\begin{array}{c}{P}_{\mathrm{issue}-2}\≈\underset{m=2}{\overset{D}{\mathrm{\Σ}}}\left(\begin{array}{c}D\\ m\end{array}\right){\left(p_0^{(L-1)}\right)}^m{(1-p_0^{(L-1)})}^{(D-m)}\\ =\underset{m=2}{\overset{D}{\mathrm{\Σ}}}\left(\begin{array}{c}D\\ m\end{array}\right){\left(\frac{N-K}{N}\right)}^{m(L-1)}{(1-{\left(\frac{N-K}{N}\right)}^{(L-1)})}^{(D-m)}\end{array}---\left(19\right)$

易于可见在K＜＜N时，可以通过为L和D选择恰当值来使两个概率相对大。例如在K＝100、N＝1000、L＝10并且D＝4时，P_issue-1＝34.87％并且P_issue-2＝50.29％。

在问题1的情况下，测量符号为零。也注意如果测量符号与一个或者多个非零源符号关联，则这些非零符号很不可能完全相互取消(见公式(16))，因此测量符号具有为非零的接近1的概率。因此如果零测量符号出现，则可以推断与它关联的所有源符号为零。

在问题2的情况下，非零源符号的两个1度测量相关如下。假设非零源符号x_i在维度d₁和d₂中具有两个1度测量。表示 $j_i^{\left(d_m\right)}=|({{\mathrm{\Π}}_d}_m\left(i\right)-1)/L|$ 并且 $l_i^{\left(d_m\right)}=\mathrm{mod}({{\mathrm{\Π}}_d}_m\left(i\right)-1,L)+1,$ 其中m＝1和2。这些两个1度测量的值具有以下关系：

$x_i=\frac{y_{j_i^{\left(d_1\right)}}^{\left(d_1\right)}}{a_{j_i^{\left(d_1\right)},l_i^{\left(d_1\right)}}^{\left(d_1\right)}}=\frac{y_{j_i^{\left(d_2\right)}}^{\left(d_2\right)}}{a_{j_i^{\left(d_2\right)},l_i^{\left(d_2\right)}}^{\left(d_2\right)}}---\left(20\right)$

显然可以从公式(20)计算x_i的值。

基于问题1和2的上述讨论，然后在步骤S42中，针对中的每个元素 ${\tilde{y}}_j^{\left(d\right)}(d=1~D,j=1~J),$ 判断是否等于0；并且如果

等于0，则设置

其中l＝1～L，其中初始化

为并且

针对稀疏信号矢量x中的每个源符号x_i，其中i＝1～N，判断是否满足以下等式：

$\frac{{\tilde{y}}_{j_i^{\left(d_m\right)}}^{\left(d_m\right)}}{a_{j_i^{\left(d_m\right)},l_i^{\left(d_m\right)}}^{\left(d_m\right)}}=\frac{{\tilde{y}}_{j_i^{\left(d_n\right)}}^{\left(d_n\right)}}{a_{j_i^{\left(d_n\right)},l_i^{\left(d_n\right)}}^{\left(d_n\right)}},{\∀d}_m=1~D,d_n=1~D,d_m\≠d_n---\left(21\right)$

如果针对一对(d_m，d_n)满足等式，则设置

然后在步骤S43中，按照以下公式通过从测量矢量y减去恢复的稀疏信号来更新

$\tilde{y}=y-A\tilde{x}---\left(22\right)$

其中

代表其中未恢复的符号设置为0的恢复的稀疏信号矢量。

这样的干扰取消操作将以两种方式帮助恢复其它符号。第一，如果测量符号让除了一个关联源符号之外的所有关联源符号恢复，则可以从测量直接估计仅有的未恢复的源符号(在干扰取消之后)。第二，如果从y减去非零源符号，则将它的关联测量的度减少1。这可能在

中引入问题1和2的新情况，因此可以再次应用上述两个规则以从重新生成更多源符号。这促使将迭代过程与干扰取消一起用来递归地恢复x。

然后在步骤S44中，将上述步骤S42和S43重复多个迭代直至满足预定条件。

预定条件可以例如包括以下任一项：

-已经达到固定迭代数目；

-成功恢复稀疏信号矢量中的所有源符号；

-

其中

和

代表在一个迭代之前和之后的

上文提出的迭代算法在这一实施例中的复杂度是O(N)。

在下文中，使用数值结果来举例说明上文提出的技术解决方案的性能。考虑以高斯随机变量为非零元素的稀疏源信号x。根据公式(14)和(15)生成测量矩阵A而高斯随机变量为非零条目。在仿真中，设置p₁＝K/N＝0.1和0.05并且维度数目D＝4。调整组长度L以获得不同概略长度M。在满足以下三个条件中的任何条件时终止迭代恢复过程：a)重新搜索T＝20个迭代；b)已经成功恢复所有源符号；以及c)

而

和为

在一个迭代之前和之后的版本。平均迭代数目根据M/N为T_ave＝6～10。

图5示出在测量数目M与Klo_g2(N/K)之间的关系。这一幅图的目的在于示出提出的技术可以实现概略长度O(Klog2(N/K))的下界。在图5中示出对于上文提出的技术解决方案，所需测量数目M随信号长度N线性地增长。这里从1000到10000设置N并且选择M以确保恢复失败概率不大于0.01。注意恢复失败出现于时，其中表示在最终迭代之后的

图5清楚地示出提出的技术解决方案可以经验地实现概略长度下界α·Klo_g2(N/K)而α约为06～0.8的很小值。

图6a和6b针对上文提出的技术解决方案示出在恢复失败概率与α(即概略长度)之间的关系。这里固定N＝4000并且从0.5～1.2改变α。为了比较，在图6中也针对基于l₁幻数和贝叶斯框架的现有恢复算法示出恢复失败概率。如可见的那样，为了实现相同恢复准确度，上文提出的技术解决方案需要比l₁幻数和贝叶斯框架少得多的测量。用线性复杂度实现这样的良好概略长度。良好经验性能和线性恢复复杂度使提出的技术成为用于利用稀疏矩阵的压缩感测的潜在和有吸引力的解决方案。

上文已经从方法观点描述本发明的技术解决方案，并且下文将从装置观点进一步描述本发明的技术解决方案。

根据本发明的一个实施例，还提供一种用于测量和恢复稀疏信号的装置。该装置包括：测量装置，用于用M×N测量矩阵A对长度为N的稀疏信号矢量x＝{x_i}执行线性测量以获得长度为M的测量矢量y，其中测量矩阵A由下式表示：

$A=\left[\begin{array}{c}{A}_1{\mathrm{\Π}}_1\\ \·\\ \·\\ \·\\ A_d{\mathrm{\Π}}_d\\ \·\\ \·\\ \·\\ A_D{\mathrm{\Π}}_D\end{array}\right],$

其中A_d∏_d是测量矩阵A的子矩阵，其中d＝1～D，∏_d是N×N随机排列矩阵，并且A_d是下式表示的J×N矩阵：

其中M＜＜N，J*L＝N，根据相同排列矩阵生成的测量符号称为一个维度，并且D是总维度数目；以及

恢复装置，用于从长度为M的测量矢量y恢复长度为N的稀疏信号矢量x＝{x_i}。

有利地，在稀疏信号是数字稀疏信号并且稀疏信号矢量x中的每个条目x_i被取自于有限集合Q＝{X₀＝0，X₁，...，X_Q-1}而X_q为非零数时，其中q＝1～Q-1，并且Q是集合的大小，恢复装置还可以包括：包括D个最大似然检测器的检测装置，用于针对D个维度分别执行最大似然检测，其中第d个最大似然检测器用于针对第d个维度执行最大似然检测；其中检测装置将上述检测重复多个迭代直至满足预定条件；以及估计装置，用于基于第D个最大似然检测器在最后迭代中的输出来估计稀疏信号矢量中的源符号。

具体而言，针对第d个维度，检测装置中的第d个最大似然检测器用于基于第d个维度的测量符号和第d个维度中的先验信息执行最大似然检测以生成稀疏信号矢量中的源符号的后验信息，其中在第d个维度中生成的后验信息用于更新第(mod(d，D)+1)个维度中的先验信息。

优选地，第d个维度中的后验信息p^(d)(x_i＝X_q)用于按照以下公式更新第(mod(d，D)+1)个维度中的先验信息

${\tilde{p}}^{(\mathrm{mod}(d,D)+1)}(x_i=X_q)=p^{\left(d\right)}(x_i=X_q)/e^{(\mathrm{mod}(d,D)+1)}(x_i=X_q),$

其中e^{(mod(d，D)+1)}(x_i＝X_q)在先前迭代中在第(mod(d，D)+1)个维度中生成，其代表在第(mod(d，D)+1)个维度中x_i为X_q的外赋概率，其中q＝0～Q-1，e^{(mod(d，D)+1)}(x_i＝X_q)在第一迭代中被初始化为1并且 $e^{(\mathrm{mod}(d,D)+1)}(x_i=X_q)=p^{(\mathrm{mod}(d,D)+1)}(x_i=X_q)/{\tilde{p}}^{(\mathrm{mod}(d,D)+1)}(x_i=X_q).$

预定条件包括以下任一项：

-已经达到固定迭代数目；

-在两个相继迭代中在第D个最大似然检测器中生成的后验概率之差在预定义阈值以下。

有利地，在稀疏信号是模拟稀疏信号时，恢复装置还可以包括：

判断装置，用于针对

中的每个元素

其中d＝1～D，j＝1～J，断

是否等于0，其中初始化

为

并且如果等于0则设置其中l＝1～L；

其中∏_d(i)是x_i在第d个排列版本中的索引，并且

是∏_d(i)的逆运算；

并且

针对稀疏信号矢量x中的每个源符号x_i，其中i＝1～N，判断是否满足以下等式：

$\frac{{\tilde{y}}_{j_i^{\left(d_m\right)}}^{\left(d_m\right)}}{a_{j_i^{\left(d_m\right)},l_i^{\left(d_m\right)}}^{\left(d_m\right)}}=\frac{{\tilde{y}}_{j_i^{\left(d_n\right)}}^{\left(d_n\right)}}{a_{j_i^{\left(d_n\right)},l_i^{\left(d_n\right)}}^{\left(d_n\right)}},{\∀d}_m=1~D,d_n=1~D,d_m\≠d_n;$

并且如果针对一对(d_m，d_n)满足该等式则设置

其中

并且 $l_i^{\left(d_m\right)}=\mathrm{mod}({{\mathrm{\Π}}_d}_m\left(i\right)-1,L)+1;$

以及

更新装置，用于按照以下公式通过从测量矢量y减去恢复的稀疏信号来更新

$\tilde{y}=y-A\tilde{x},$

其中代表其中未恢复的符号设置为0的恢复的稀疏信号矢量；

其中判断装置将上述判断重复多个迭代直至满足预定条件。

预定条件可以例如包括以下任一项：

-已经达到固定迭代数目；

-成功恢复稀疏信号矢量中的所有源符号；

-

其中

和代表在一个迭代之前和之后的

应当注意给出上文描述的实施例用于描述而不是限制实施例，并且将理解可以如本领域技术人员容易理解的那样求助于修改和变化而未脱离本发明的精神实质和范围。这样的修改和变化视为在本发明和所附权利要求的范围内。本发明的保护范围由所附权利要求限定。此外，不应解释权利要求中的标号为对权利要求的限制。使用动词“包括”及其变形未排除除了在权利要求中声明的单元或者步骤之外的单元或者步骤的存在。在单元或者步骤之前的不定冠词“一个”未排除多个这样的单元或者步骤的存在。

Claims (14)

1.一种用于测量和恢复稀疏信号的方法，包括：

a.用M×N测量矩阵A对长度为N的稀疏信号矢量x＝{x_i}执行线性测量以获得长度为M的测量矢量y，其中所述测量矩阵A由下式表示：

$A=\left[\begin{array}{c}{A}_1{\mathrm{\Π}}_1\\ \·\\ \·\\ \·\\ A_d{\mathrm{\Π}}_d\\ \·\\ \·\\ \·\\ A_D{\mathrm{\Π}}_D\end{array}\right],$

其中A_d∏_d是所述测量矩阵A的子矩阵，其中d＝1～D，∏_d是N×N随机排列矩阵，并且A_d是下式表示的J×N矩阵：

其中M＜＜N，J*L＝N，根据相同排列矩阵生成的测量符号称为一个维度，并且所述D是总维度数目；并且

b.从所述长度为M的测量矢量y恢复所述长度为N的稀疏信号矢量x＝{x_i}。

2.根据权利要求1所述的方法，其中所述稀疏信号是数字稀疏信号，并且所述稀疏信号矢量x中的每个条目x_i被取自于有限集合Q＝{X₀＝0，X₁，...，X_Q-1}而X_q为非零数，其中q＝1～Q-1，并且Q是所述集合的大小，所述步骤b包括：

b1.使用D个最大似然检测器以针对D个维度分别执行最大似然检测，其中第d个最大似然检测器用于针对第d个维度执行最大似然检测；

b2.将步骤b1重复多个迭代直至满足预定条件；

b3.基于第D个最大似然检测器在最后迭代中的输出来估计所述稀疏信号矢量中的源符号。

3.根据权利要求2所述的方法，其中所述步骤b1包括：

-针对所述第d个维度，基于所述第d个维度的所述测量符号和所述第d个维度中的先验信息执行最大似然检测以生成所述稀疏信号矢量中的所述源符号的后验信息，其中在所述第d个维度中生成的所述后验信息用于更新第(mod(d，D)+1)个维度中的先验信息。

4.根据权利要求3所述的方法，其中所述第d个维度中的所述后验信息p^(d)(x_i＝X_q)用于按照以下公式更新所述第(mod(d，D)+1)个维度中的所述先验信息 ${\tilde{p}}^{(\mathrm{mod}(d,D)+1)}(x_i=X_q):$

${\tilde{p}}^{(\mathrm{mod}(d,D)+1)}(x_i=X_q)=p^{\left(d\right)}(x_i=X_q)/e^{(\mathrm{mod}(d,D)+1)}(x_i=X_q),$

其中所述e^{(mod(d，D)+1)}(x_i＝X_q)在先前迭代中在所述第(mod(d，D)+1)个维度中生成，其代表在所述第(mod(d，D)+1)个维度中x_i为X_q的外赋概率，其中q＝0～Q-1，所述e^{(mod(d，D)+1)}(x_i＝X_q)在第一迭代中被初始化为1并且在后继迭代中被更新为 $e^{(\mathrm{mod}(d,D)+1)}(x_i=X_q)=p^{(\mathrm{mod}(d,D)+1)}(x_i=X_q)/{\tilde{p}}^{(\mathrm{mod}(d,D)+1)}(x_i=X_q).$

5.根据权利要求1所述的方法，其中在所述稀疏信号是模拟稀疏信号时，所述方法b包括：

i.针对

中的每个元素

其中d＝1～D，j＝1～J，判断所述

是否等于0，其中初始化为

并且

是

的第d个维度中的第j个元素；

ii.如果所述

等于0，则设置

其中l＝1～L；

其中∏_d(i)是x_i在第d个排列版本中的索引，并且

是∏_d(i)的逆运算；

并且所述方法还包括：

u.针对所述稀疏信号矢量x中的每个源符号x_i，其中i＝1～N，判断是否满足以下等式：

$\frac{{\tilde{y}}_{j_i^{\left(d_m\right)}}^{\left(d_m\right)}}{a_{j_i^{\left(d_m\right)},l_i^{\left(d_m\right)}}^{\left(d_m\right)}}=\frac{{\tilde{y}}_{j_i^{\left(d_n\right)}}^{\left(d_n\right)}}{a_{j_i^{\left(d_n\right)},l_i^{\left(d_n\right)}}^{\left(d_n\right)}},\∀d_m=1~D,d_n=1~D,d_m\≠d_n;$

v.如果针对一对(d_m，d_n)满足所述等式，则设置

其中 $j_i^{\left(d_m\right)}=|({{\mathrm{\Π}}_d}_m\left(i\right)-1)/L|$ 并且 $l_i^{\left(d_m\right)}=\mathrm{mod}({{\mathrm{\Π}}_d}_m\left(i\right)-1,L)+1;$

并且在上述步骤之后，所述方法还包括：

p.按照以下公式通过从所述测量矢量y减去所述恢复的稀疏信号来更新

$\tilde{y}=y-A\tilde{x},$

其中

代表其中未恢复的符号设置为0的所述恢复的稀疏信号矢量；

q.将上述步骤重复多个迭代直至满足预定条件。

6.根据权利要求2所述的方法，其中所述预定条件包括以下任一项：

-已经达到固定迭代数目；

-在两个相继迭代中在第D个最大似然检测器中生成的后验概率之差在预定义阈值以下。

7.根据权利要求5所述的方法，其中所述预定条件包括以下任一项：

-已经达到固定迭代数目；

-成功恢复所述稀疏信号矢量中的所有源符号；

-

其中

和

代表在一个迭代之前和之后的

8.一种用于测量和恢复稀疏信号的装置，包括：

测量装置，用于用M×N测量矩阵A对长度为N的稀疏信号矢量x＝{x_i}执行线性测量以获得长度为M的测量矢量y，其中所述测量矩阵A由下式表示：

$A=\left[\begin{array}{c}{A}_1{\mathrm{\Π}}_1\\ \·\\ \·\\ \·\\ A_d{\mathrm{\Π}}_d\\ \·\\ \·\\ \·\\ A_D{\mathrm{\Π}}_D\end{array}\right],$

其中A_d∏_d是所述测量矩阵A的子矩阵，其中d＝1～D，∏_d是N×N随机排列矩阵，并且A_d是下式表示的J×N矩阵：

其中M＜＜N，J*L＝N，根据相同排列矩阵生成的测量符号称为一个维度，并且所述D是总维度数目；以及

恢复装置，用于从所述长度为M的测量矢量y恢复所述长度为N的稀疏信号矢量x＝{x_i}。

9.根据权利要求8所述的装置，其中所述稀疏信号是数字稀疏信号，并且所述稀疏信号矢量x中的每个条目x_i被取自于有限集合Q＝{X₀＝0，X₁，...，X_Q-1}而X_q为非零数，其中q＝1～Q-1，并且Q是所述集合的大小，所述恢复装置包括：

包括D个最大似然检测器的检测装置，用于针对D个维度分别执行最大似然检测，其中第d个最大似然检测器用于针对第d个维度执行最大似然检测；

其中所述检测装置将上述检测重复多个迭代直至满足预定条件；

估计装置，用于基于第D个最大似然检测器在最后迭代中的输出来估计所述稀疏信号矢量中的源符号。

10.根据权利要求9所述的装置，其中针对所述第d个维度，所述检测装置中的所述第d个最大似然检测器用于基于所述第d个维度的所述测量符号和所述第d个维度中的先验信息执行最大似然检测以生成所述稀疏信号矢量中的所述源符号的后验信息，其中在所述第d个维度中生成的所述后验信息用于更新第(mod(d，D)+1)个维度中的先验信息。

11.根据权利要求10所述的装置，其中所述第d个维度中的所述后验信息p^(d)(x_i＝X_q)用于按照以下公式更新所述第(mod(d，D)+1)个维度中的所述先验信息 ${\tilde{p}}^{(\mathrm{mod}(d,D)+1)}(x_i=X_q):$

${\tilde{p}}^{(\mathrm{mod}(d,D)+1)}(x_i=X_q)=p^{\left(d\right)}(x_i=X_q)/e^{(\mathrm{mod}(d,D)+1)}(x_i=X_q),$

其中所述e^{(mod(d，D)+1)}(x_i＝X_q)在先前迭代中在所述第(mod(d，D)+1)个维度中生成，其代表在所述第(mod(d，D)+1)个维度中x_i为X_q的外赋概率，其中q＝0～Q-1，所述e^{(mod(d，D)+1)}(x_i＝X_q)在第一迭代中被初始化为1并且在后继迭代中被更新为 $e^{(\mathrm{mod}(d,D)+1)}(x_i=X_q)=p^{(\mathrm{mod}(d,D)+1)}(x_i=X_q)/{\tilde{p}}^{(\mathrm{mod}(d,D)+1)}(x_i=X_q).$

12.根据权利要求8所述的装置，其中在所述稀疏信号是模拟稀疏信号时，所述恢复装置包括：

判断装置，用于针对

中的每个元素

其中d＝1～D，j＝1～J，判断所述是否等于0，其中初始化为

并且是

的第d个维度中的第j个元素；

并且如果所述

等于0则设置其中l＝1～L；

其中∏_d(i)是x_i在第d个排列版本中的索引，并且

是∏_d(i)的逆运算；

并且

针对所述稀疏信号矢量x中的每个源符号x_i，其中i＝1～N，判断是否满足以下等式：

$\frac{{\tilde{y}}_{j_i^{\left(d_m\right)}}^{\left(d_m\right)}}{a_{j_i^{\left(d_m\right)},l_i^{\left(d_m\right)}}^{\left(d_m\right)}}=\frac{{\tilde{y}}_{j_i^{\left(d_n\right)}}^{\left(d_n\right)}}{a_{j_i^{\left(d_n\right)},l_i^{\left(d_n\right)}}^{\left(d_n\right)}},{\∀d}_m=1~D,d_n=1~D,d_m\≠d_n;$

并且如果针对一对(d_m，d_n)满足所述等式则设置

其中并且 $l_i^{\left(d_m\right)}=\mathrm{mod}({{\mathrm{\Π}}_d}_m\left(i\right)-1,L)+1;$

更新装置，用于按照以下公式通过从所述测量矢量y减去所述恢复的稀疏信号来更新

$\tilde{y}=y-A\tilde{x},$

其中

代表其中未恢复的符号设置为0的所述恢复的稀疏信号矢量；

其中所述判断装置将上述判断重复多个迭代直至满足预定条件。

13.根据权利要求9所述的装置，其中所述预定条件包括以下任一项：

-已经达到固定迭代数目；

-在两个相继迭代中在第D个最大似然检测器中生成的后验概率之差在预定义阈值以下。

14.根据权利要求12所述的装置，其中所述预定条件包括以下任一项：

-已经达到固定迭代数目；

-成功恢复所述稀疏信号矢量中的所有源符号；

-

其中

和

代表在一个迭代之前和之后的

Images