压缩感知——简介
 
压缩感知，compressed sensing又称compressed sampling，是在采样过程中完成了数据压缩的过程。
 压缩感知在信号采样的过程中，用很少的采样点，实现了和全采样一样的效果。
 
压缩感知——信号采样
 
学过通信原理或信号与系统的都知道奈奎斯特采样定理，即想让采样之后的数字信号完整保留原始信号中的信息，采样频率必须大于信号中最高频率的2倍。原因是时域以τ为间隔进行采样，频域会以1/τ为周期发生周期延拓。那么如果采样频率低于两倍的信号最高频率，信号在频域频谱搬移后就会发生混叠。
 2004年，Candes，陶哲轩和Donoho提出了压缩感知理论，该理论认为：如果信号是稀疏的，那么它可以由远低于采样定理要求的采样点重建恢复。
 突破的关键就在于采样的方式。当我们说“采样频率”的时候，意味着做的是等间距采样，数字信号领域通常都是做等间距采样，也服从奈奎斯特采样定律。
 
 但是压缩感知采用不等间距采样，比如采用随机亚采样（图b上方的红点），那么这时候频域就不再是以固定周期进行延拓了，而是会产生大量不相关（incoherent）的干扰值。如图c，最大的几个峰值还依稀可见，只是一定程度上被干扰值覆盖。这些干扰值看上去非常像随机噪声，但实际上是由于三个原始信号的非零值发生能量泄露导致的（不同颜色的干扰值表示它们分别是由于对应颜色的原始信号的非零值泄露导致的）
 之所以随机亚采样会有这样的效果，可以理解成随机采样使得频谱不再是整齐地搬移，而是一小部分一小部分胡乱地搬移，频率泄露均匀地分布在整个频域，因而泄漏值都比较小，从而有了恢复的可能。
 
压缩感知——信号恢复
 
压缩感知——两个前提条件
 
刚刚的例子之所以能够实现最终信号的恢复，是因为它满足了两个前提条件：
 
这个信号在频域只有3个非零值，所以可以较轻松地恢复出它们。
采用了随机亚采样机制，因而使频率泄露均匀地分布在整个频域。
 这两点对应了CS的两个前提条件——稀疏性（sparsity）、不相关性（incoherence）。
 稀疏性可以简单直观地理解为：若信号在某个域中只有少量非零值，即信号在某个域中非零点远远小于信号总点数，那么它在该域稀疏，该域也被称为信号的稀疏域。
 因此，第一个前提条件要求信号必须在某一个变换域具有稀疏性。比如例子中，信号在频域是稀疏的，因而可以通过所述的重建方法轻松地在稀疏域（频域）复原出原信号。
 然而通常信号在变换域中不会呈现完全的稀疏性。其实只要它近似满足稀疏性，即大部分值趋于零，只有少量大的非零值，就可以认为它是可压缩信号，可以对它进行CS亚采样。
 对于之前讲的例子，如果它在频域中不稀疏，我们可以做DWT、DCT等，找到它的稀疏变换。
 
前提条件1——稀疏性
 
前提条件2——非相关性/等距约束性
 
在讲第二个前提条件之前，需要引入必要的数学表达。
 
 上面这张图把亚采样的过程用矩阵的方式表达出来：
 如图，x是为长度N的一维信号，也就是原信号，稀疏度为k。此刻它是未知的。
 Φ为观测矩阵，对应着亚采样这一过程。它将高维信号x投影到低维空间，是已知的。
 y=Φx为长度M的一维测量值，也就是亚采样后的结果。显然它也是已知的。
 因此，压缩感知问题就是在已知测量值y和测量矩阵Φ的基础上，求解欠定方程组y=Φx得到原信号x。
 然而，一般的自然信号x本身并不是稀疏的，需要在某种稀疏基上进行稀疏表示。令x=Ψs，Ψ为稀疏基矩阵，s为稀疏系数。
 于是最终方程就变成了：y=ΦΨs。已知y、Φ、Ψ，求解s。
 y=ΦΨs有点长，我们把ΦΨ合并成一个矩阵，称之为传感矩阵。即令Θ=ΦΨ ，则y=ΘS。问题即为，已知y和Θ，求解S。
 对应到一开始的例子中：
 x就是三个正弦信号叠加在一起的原信号a；
 稀疏矩阵Ψ就是傅里叶变换，将信号变换到频域S；
 观测矩阵Φ就是我们采用的随机亚采样方式；
 y就是最终的随机亚采样的结果c。
 
 求解出S后，由x=Ψs即可得到恢复出的原信号x。
 然而在正常情况下，方程的个数远小于未知数的个数，方程是没有确定解的，无法重构信号。但是，由于信号是K稀疏，如果上式中的Φ满足有限等距性质(RIP)，则K个系数就能够从M个测量值准确重构（得到一个最优解）。
 陶哲轩和Candès证明了RIP是观测矩阵要满足的准确要求。但是，要确认一个矩阵是否满足RIP非常复杂。于是Baraniuk证明：RIP的等价条件是观测矩阵和稀疏表示基不相关（incoherent）。即压缩感知的第二个前提条件。
 
 
 那怎样找到不相关的观测矩阵呢？陶哲轩和Candès又证明： 独立同分布的高斯随机测量矩阵可以成为普适的压缩感知测量矩阵。
 于是满足高斯分布的随机测量矩阵就成了CS最常用的观测矩阵。
 对于二维信号，往往就采用二维高斯随机测量采样矩阵对图像进行亚采样。
 对于一维信号，采用前文提到的随机不等间距的亚采样即可。
 
压缩感知——重建方法
 
信号采样后需要将其恢复。CS的重建也就是求解欠定方程组y=ΘS的方法。这是一个零范数（l0）最小化问题，是一个NP完全问题（没有快速解法的问题），因此往往转换成一范数（l1）最小化的求解，或者用一些近似估计的算法。
 匹配追踪是一种典型的算法。以上文中的例子为例：
 (1) 由于原信号的频率非零值在亚采样后的频域中依然保留较大的值，其中较大的两个可以通过设置阈值，检测出来（图a）。
 (2) 然后，假设信号只存在这两个非零值（图b），则可以计算出由这两个非零值引起的干扰（图c）。
 (3) 用a减去c，即可得到仅由蓝色非零值和由它导致的干扰值(图d），再设置阈值即可检测出它，得到最终复原频域（图e）
 (4) 如果原信号频域中有更多的非零值，则可通过迭代将其一一解出。
 
 
扩展：图像压缩与压缩感知
 
信号的稀疏性已经在图像压缩领域有了很广泛的应用。利用信号的稀疏性，可以对信号进行压缩。如图像压缩领域的JPEG格式，就是将图像变换到离散余弦域，得到近似稀疏矩阵，只保留较大的值，从而实现压缩。
 图像压缩和压缩感知这两个概念有着本质上的区别。
 图像压缩是先进行了全采样，然后再变换域丢弃小系数，完成压缩；
 而压缩感知不同，它的思想其实从图像压缩中借鉴了很多：既然全采样了还要再丢弃，我们为什么不能直接少采样一些点？因此，压缩感知直接进行了亚采样，然后再用算法消除亚采样导致的伪影。可以说，压缩感知直接在采样时就完成了压缩。
 
 这种直接减少采样点，采集完后重建的方式，相比于图像压缩的全采用之后再压缩的方式，对于很多情形，比如照相机拍摄照片，并没有优势。因为所有像素的采集在一瞬间就都完成了。但是对于一些采集比较慢的情形，比如核磁共振成像，CS就可以发挥巨大优势。原本一副MRI图像常常需要几十秒，速度慢也是MRI的一大缺陷。而应用CS技术后，只需要采集全采样几分之一的数据，就可以重建出原图。这样就可以把成像速度提高好几倍，同时对图像质量影响不大。
 
总结
 
什么是压缩感知：
 如果一个信号在某个变换域是稀疏的，那么就可以用一个与变换基不相关的观测矩阵将变换所得高维信号投影到一个低维空间上，然后通过求解一个优化问题就可以从这些少量的投影中以高概率重构出原信号。
 压缩感知理论主要包括三部分：
 （1）信号的稀疏表示；
 （2）设计测量矩阵，要在降低维数的同时保证原始信号x的信息损失最小；
 （3）设计信号恢复算法，利用M个观测值无失真地恢复出长度为N的原始信号。
 
从06年提出至今，已经发展出了很多算法，原来的基于l1 minimization的BP算法很慢，现在都是快速算法，而且求解算法也从纯优化方面扩展到了estimation方面，有很多基于贝叶斯估计的方法了，目前最火的也是Donoho他们组搞得AMP算法，是用Graph model里面的message passing算法通过近似求解MMSE（MAP）解。在测量矩阵方面，也已经设计出了各种矩阵，除了i.i.d. Gaussian的矩阵还有很多正交的矩阵，比如partial random DFT/DCT 矩阵。对信号的要求也从稀疏变成了存在某种结构，比如low rank，group sparse等等。（2017年进展）
 
压缩感知和矩阵填充
 
矩阵填充的要领是通过低秩矩阵中的已知要素还原出该矩阵的其他未知要素的进程。这几年，关于矩阵填充方法的理论研究成为压缩感知技术的一个研究热点。在实际的应用领域中涉及对高维数据的分析与处理，可以运用矩阵填充的方法来解决。其过程主要是：通过观测到的局部数据来准确填充缺失数据,从而获得完整数据矩阵的过程。
 压缩感知和矩阵填充都是稀疏约束下的反问题，压缩感知利用信号本身或在变换域中的稀疏约束性求解欠定方程，矩阵填充利用矩阵的低秩约束性求解欠定方程。
 压缩感知理论的核心问题是信号的稀疏表示、观测矩阵的设计和重构算法，信号本身或在变换域中的系数越稀疏，观测矩阵和稀疏基构成的压缩感知矩阵的受限等距常数越小，则压缩感知的性能越好。
 矩阵填充理论的核心问题是矩阵的低秩特性、非相干特性和重构算法，寻找性能良好的重构算法一直是矩阵填充理论中的一个研究重点。此外，压缩感知的应用领域已经拓展得较为广泛，但矩阵填充的应用尚处于起步阶段，挖掘矩阵填充的应用，进而将矩阵填充和压缩感知结合起来进行应用方面的探索，是非常重要和有意义的课题。
 压缩感知的性能取决于3个要素：信号的稀疏性、压缩感知矩阵的非相干性和重构算法的快速有效性。
 矩阵填充性能也取决于3个要素：矩阵的低秩性、矩阵的不相关性和重构算法的快速有效性。
 
低秩矩阵
 
还记得我们怎么手工求矩阵的秩吗？为了求矩阵A的秩，我们是通过矩阵初等变换把A化为阶梯型矩阵，若该阶梯型矩阵有r个非零行，那A的秩rank(A)就等于r。从物理意义上讲，矩阵的秩度量的就是矩阵的行列之间的相关性。如果矩阵的各行或列是线性无关的，矩阵就是满秩的，也就是秩等于行数。回到上面线性方程组来说吧，因为线性方程组可以用矩阵描述嘛。秩就表示了有多少个有用的方程了。上面的方程组有3个方程，实际上只有2个是有用的，一个是多余的，所以对应的矩阵的秩就是2了。
 OK。既然秩可以度量相关性，而矩阵的相关性实际上就表示了矩阵的结构信息。如果矩阵之间各行的相关性很强，那么就表示这个矩阵实际可以投影到更低维的线性子空间，也就是用几个向量就可以完全表达了，它就是低秩的。所以我们总结的一点就是：如果矩阵表达的是结构性信息，例如图像、用户-商品推荐表等等，那么这个矩阵各行之间存在这一定的相关性，那这个矩阵一般就是低秩的。
 如果X是一个m行n列的数值矩阵，rank(X)是X的秩，假如rank (X)远小于m和n，则我们称X是低秩矩阵。低秩矩阵每行或每列都可以用其他的行或列线性表出，可见它包含大量的冗余信息。利用这种冗余信息，可以对缺失数据进行恢复，也可以对数据进行特征提取。
 
应用
 
矩阵填充越来越多的应用在协同滤波、系统识别、无线传感网络、视频去噪、人脸识别、医学成像等领域，正在发挥着巨大的作用。特别是在室内定位中的 应用越来越多，是当下的研究热点之一。
 
参考网址：
 形象易懂讲解算法II——压缩感知（非常好）
 AndyJee：浅谈压缩感知系列——博客园（非常丰富）
 压缩感知理论
 机器学习——低秩矩阵分解中低秩的意义、矩阵填补、交叉验证

压缩感知理论及其算法研究报告
 
摘 要:随着信息技术的发展，人们对信息的巨量需求以及硬件的发展缓慢造成了信号采样，传输和存储的巨大压力。如何解决在现有的硬件基础上传输大量的信息成为热点研究的内容。近年来压缩感知的出现为缓解这些压力提供了解决的办法。本文综述了压缩感知的理论框架及关键的技术问题，并着重介绍了压缩感知稀疏重构中的主流贪婪算法，通过算法实验分析了各种算法的重构性能。
 关键词：压缩感知 贪婪算法 稀疏重构
 1.引言
 传统的信号采样定律-那奎斯特采样定律定理：为了不失真地恢复模拟信号，采样频率应该不小于模拟信号频谱中最高频率的2倍[1]。这样对于系统处理信息的硬件需求提出了很高的要求，同时在实际应用当中为了节约存储空间和降低传输成本，需要对采集的数据进行压缩处理，这样会造成大量采集的数据浪费。因而压缩感知技术应运而生，打破了传统的信号采样定理，从不同的角度解决了信号采样的问题。使得在保证信息不损失的情况下,用远低于奈奎斯特采样定理要求的速率采样信号。压缩感知理论指出只要信号是可压缩的或在某个变换域是稀疏的,那么就可以用一个与变换基不相关的观测矩阵将变换所得高维信号投影到一个低维空间上,然后通过求解一个优化问题就可以从这些少量的投影中以高概率重构出原信号[2]。在该理论框架下,采样速率不决定于信号的带宽,而决定于信息在信号中的结构和内容。压缩感知理论使得采样和计算的成本大大降低。
 当前主流的压缩感知重构算法主要包括三类：凸优化方法,贪婪算法和基于贝叶斯框架提出的算法[3]。本文主要以压缩感知重构算法为主线，介绍了主流的贪婪追踪类算法包括正交匹配追踪(OMP)算法，正则化正交匹配追踪(ROMP)算法，分段正交匹配追踪(STOMP)算法，稀疏度自适应匹配追踪(SAMP)算法，并且对这些算法的优缺点进行了比较通过仿真实验分析了各种算法的重构性能。
 2.压缩感知基本理论
 压缩感知理论主要包括信号的稀疏表示，编码测量和信号重构算法三个方面。信号的稀疏表示是压缩感知的先验条件，信号的稀疏表示是将信号投影到正交变换基时，绝大多数的变换稀疏的绝对值很小，所得到的变换向量是稀疏的或者是近似稀疏的。任意的N维信号都可以通过某个稀疏矩阵线性表示。例如x为N维的信号，Ψ是x对应的稀疏矩阵，则x可以表示为：x = ∑_(i=1)^N▒〖θ_i ψ_i 〗其中，Ψ是N列ψ_i组成的矩阵，θ_i是x在ψ_i下的投影系数，θ是投影系数向量，θ=Ψ^Tx。稀疏矩阵一般根据信号本身特点灵活选取，常见的是离散余弦变换基，快速傅里叶变换基，离散小波变换基。在编码测量中首先选定一个平稳的，与稀疏基Ψ不相关的M × N 维的观测矩阵Ф，对θ进行观测得到观测集合Y = Фx = ФΨθ。令A=ФΨ为M × N 的矩阵，称为感知矩阵。Y可以看作是稀疏信号θ关于测量矩阵A的测量值,上式整体可以表示为图(1)。
 
 
在压缩感知的整个过程中,测量矩阵的设计是一个关键步骤。测量矩阵性质的好坏,关系到能否达到压缩的目的,同时又直接关系到信号能否被精确重构。设计一个合适的观测矩阵应该既能达到压缩采样的目的,同时又可以保证信号可以无失真的重构。有限等距性质[4]在理论上较好的解决了这个问题，只要感知矩阵A能够满足RIP条件,那么信号可以由少量的测量值经过重构算法精确的恢复出来,也就是说,理论上我们可以设计一个测量矩阵使得感知矩阵A满足RIP规则，这样既可以达到压缩采样的目的,又能保证信号无失真的恢复出来。RIP规则的数学表达描述为：设A=ФΨ为M × N 的矩阵，假设一个常数δ_k,使得对于任意向量s和所有的矩阵A_k,满足以下关系4
 (1-δ_k)〖||s||〗_2≤〖||A_k s||〗_2≤(1+δ_k)〖||s||〗2 (2-1)
 其中A_k是A子矩阵,大小为M×K,有限等距常数为δ_k∈(0,1)。
 如果A满足约束等距原则，保证了信号恢复的唯一性。实际上要直接验证矩阵是否满足RIP条件是一件很困难的事情。在实际应用中,我们可以用RIP准则的一种等价情况,即非相干性来指导测量矩阵的设计。非相干性指测量矩阵中的行向量不能被稀疏矩阵线性表出同理稀疏矩阵中的列向量也不能被测量矩阵中的任意行向量线性表出。相干性的度量由相干度[5]如图(2-2)所示，关系数旳取值范围为U∈(1，√N)
 U(Ф,Ψ)=√Nmax{|Ф_k,Ψ_J|} (2-2)
 Donoho等人在文献[6]中指出服从高斯分布的随机矩阵可以高概率满足不相关性,对于一个大小为M×N的随机高斯矩阵Ф,Ф中每个值满足均值为0，方差为1/M的高斯分布，即Ф(i，j) ~ N(0,1/M)。可以证明当M>cKlog(N/K)时，A = ФΨ在很大概率下能满足RIP条件。而且随机高斯矩阵与大多数固定正交基构成的矩阵不相关,因此随机高斯测量矩阵满足理论上的最优性。目前大多数情况下都采用随机高斯矩阵作为压缩感知的测量矩阵，本文实验所用到的观测矩阵也是随机高斯矩阵。当选取好观测矩阵，压缩感知问题转化为求解(2-1)式的最优l_0范数问题
 min〖||θ||〗_0 s.t. Aθ = Y (2-3)
 如果得到x的稀疏表示θ，可以进一步由变换基Ψ通过下式(2-2)重构原始信号
 x = Ψθ (2-4)
 由于矩阵A的维度为M × N(M << N)，所以方程(1)有无穷多解，通过贪婪算法可以逐步逼近最优解，直到求出原始信号。最早提出的是匹配追踪算法(MP),MP算法的基本思想是在每一次迭代过程中。从感知矩阵中选择与信号最匹配的原子来进行稀疏逼近求出余量，在稀疏度已知的情况下继续迭代选出与余量最匹配的原子。最匹配是指当前余量与原子的内积最大。经过数次迭代，该信号便可由这些原子线性表示，但是当前的余量仅与当前的原子正交而不是与已选定的所有的原子正交使，得每次迭代的结果可能是次最优的往往需要迭代多次。下面介绍四种主流的贪婪类重构算法并分析它们的优缺点进行比较。
 3.正交匹配追踪类算法
 3.1正交匹配追踪(OMP)算法
 OMP算法作为MP算法的延申，仍然沿用了MP算法中原子选择的标准，不同的是OMP算法利用Gram-Schmidt正交化对已选定的原子进行正交化处理，再将信号在这些正交原子构成的张量空间投影，得到信号在选定原子上的分量和余量，然后用相同的方法迭代分解余量。通过每次对所选原子的正交化处理保证了迭代的最优性，从而减少了迭代的次数[7]。
 OMP的具体步骤如下：
 (1)令初始余量r_0 = Y,稀疏度为K，索引值集合J为空集，支撑集合Λ为空集；
 (2)计算相关系数u(余量与原子的内积)，并将u中最大值对应的索引值存入J；
 (3)更新支撑集合 ，将更新索引对应的原子存入Λ；
 (4)利用最小二乘法得到重建信号，同时对余量进行跟新；
 (5)若迭代次数小于K，r = r_new,n = n+1,转到第二步继续迭代；否则，停止迭代。
 为了说明OMP算法的重建性能，利用MATLAB R2016a作为平台，分两次实验验证。第一次实验，假定信号是稀疏的并且稀疏度为10，定常信号256 × 1的列向量中随机选取10个随机数随机排列到定长信号中。观测矩阵为128×256高斯随机矩阵。初始信号和恢复信号之间的误差用二者相减的二范数表示。实验结果如图(2)所示，实验误差为o(10e-15 )
 
 
 
 
第二次实验，测试对象为256×256lena图像。观测矩阵采用随机高斯矩阵，M和N表示观测矩阵的行数和列数，M/N表示压缩比0.5，实验结果如图(3)，PSRN值为26.5536。
 

 
 
OMP算法虽然保证了每次迭代的最优性，减少了迭代的次数。但是，它每次迭代中仅选取一个原子来跟新原子的集合，这样必然会付出巨大的重建时间代价。OMP算法首次把最小二乘法引入压缩感知重建中,用正交化的思想来计算重建信号使得结果更加准确,这是压缩感知重建算法取得重大研究进展的一个标志[8]。
 3.2正则化正交匹配追踪(ROMP)算法
 ROMP算法首先根据相关原则进行原子的一次筛选，通过求余量r与测量矩阵Ф中各个原子之间的内积的绝对值，来计算相关系数u，并按照此方法筛选出的K个原子的索引值存到候选集J中以便进行原子的二次筛选。
 ROMP算法采用正则化过程进行原子的二次筛选，将J中索引值对应的原子的相关系数分成若干，要求分组的原子在各自所在子集内的原子同误差向量的内积的最大值与最小值的比值在两倍以内[6]。数学表达式为(3-1)
 |u(i)| ≤ 2|u(j)|, i,j ∈J (3-1)
 然后选择能量最大的一组相关系数对应的原子索引值存入J_0中，该正则化过程可以使得ROMP算法最多经过K次迭代便可得到一个原子数|Λ|小于2K的支撑集Ф_Λ用于重建信号，对于没有选入支撑集的原子，正则化过程则能保证它们的能量一定远小于被选原子的能量，是一种简单有效的原子筛选方法。
 ROMP的步骤：
 (1)初始化余量r_0 = Y，估计信号稀疏度为K，迭代次数n = 1，索引值Λ为空集，J为空集；
 (2)计算相关系数u，并从u中寻找K个最大值对应的索引值存入J中；
 (3)对J中索引值对应原子的相关系数进行正则化，并将正则化的结果存入J_0；
 (4)更新支撑集Ф_Λ，其中Λ=Λ∪J_0；
 (5)利用最小二乘法得到重建信号，同时对余量进行跟新；
 (6)若|Λ| 2K，则停止迭代，否则令r =r_new，n = n+1，转到步骤(2)继续迭代。
 为了说明ROMP算法的重建性能，分两次实验实现算法重构，利用MATLAB R2016a作为平台。实验1假定信号是稀疏的并且稀疏度为10，定常信号256 × 1的列向量中随机选取10个随机数随机排列到定长信号中。观测矩阵为高斯随机矩阵。初始信号和恢复信号之间的误差用二者相减的二范数表示。实验结果如图(4)误差为o(10 )
 
实验2假定信号是四个余弦函数式表示
 x=0.3cos(2π50t)+0.6cos(2π100t)+0.1cos(2π200t)+0.9cos(2π400t)
 我们运用快速傅里叶变换对一维信号 x 进行变换，稀疏度为7，测量矩阵为64×256高斯矩阵实验结果如图(5)，误差为o(10 )。
 
ROMP算法在OMP算法的基础上加入了正则化方法,以实现一次迭代选择多个原子的目的,从而提高重建速度，减少了算法的复杂度。但其代价是需要预估计信号的稀疏度,并且需要较多的采样数据
 3.3分段正交匹配追踪(STOMP)算法
 STOMP 算法采用分阶段的思想首先根据相关原则来筛选原子，利用阈值的方法从原子集合中选择和迭代余量匹配的原子，与OMP 算法不同的是，它并不是每次固定选择一个匹配原子，而是给定标准为大于门限值t_s δ_t的原子。δ_t为规范噪音水平, δ_t=〖||r_t ||〗2/√M。r_t是上一次跟新的余量值[9]。利用此标准可以一次找到多个原子，减少了匹配的次数，提高了追踪的效率; 然后更新支撑集和原子，并用最小二乘法求得近似解，同时完成对余量的更新。
 STOMP算法步骤：
 (1)初始化余量r_0 = Y，迭代次数默认为10，门限参数t_s默认为2.5。计数器t=1；
 (2)计算感知矩阵各列原子与余量的内积，选择大于门限值t_s δ_t的感知矩阵对应的原子列向量存入集合J_0 ；
 (3)跟新支撑集Λ_t = Λ(t-1) ∪ J_0，在支撑集上利用最小二乘法求出跟新余量；
 (4)若t值小于迭代次数s，返回(2)继续执行，否则退出循环。
 为了说明STOM算法的重建性能，利用MATLAB R2016a作为平台。假定信号是稀疏的并且稀疏度为10，定常信号256 × 1的列向量中随机选取10个随机数随机排列到定长信号中。观测矩阵为高斯随机矩阵。初始信号和恢复信号之间的误差用二者相减的二范数表示。实验结果如图(6)所示
 
 
 
 
STOMP 算法将 OMP 算法进行了一定程度的简化，提高了计算速度，但是由于其在每次迭代的过程中寻找的都不是信号的最佳表示，致使重构的精度降低，导致在实际中其重构的信号的精度远不如OMP 算法重构的信号的精度。门限参数和默认迭代次数的设置很大程度决定了算法的精度，使得算法的灵活度差。
 3.4稀疏度自适应匹配追踪(SAMP)算法
 以上算法均建立在稀疏度已知的情况下才能有效,但在实际情况中,信号的稀疏度信息往往事先是不可得的,SAMP相比较其他算法最吸引的就是可以不需要知道稀疏度的先验信息 ,这使得SAMP重建算法在实际应用中更加广泛[10]。
 SAMP算法中采用转换阶段(stage) 的方式逐步增加该原子数，将同一个迭代过程分成多个阶段，设置一个可变步长(size)代替所选原子数目，相邻两个阶段所对应的支撑集的大小之差即为当前步长，随着步长的增加和支撑集的不断增大，实现了在未知稀疏度的前提下步长逐步逼近稀疏度K进而实现精确重建出原始信号的目的[11]。SAMP算法引入了回溯的思想，每次迭代都重新评估原子的有效性。
 SAMP算法步骤：
 (1)初始化余量r_0= Y，支撑集F_0为空集，候选集S_0为空集，初始步长为s，步长增量初始倍数j=1；
 (2)选择感知矩阵内的原子和余量内积最大的s个列向量，将对应的原子列向量放入候选集s_k中；
 (3)合并支撑集F_k和候选集s_k得到跟新后的支撑集C_k；
 (4)以C_k的原子为基准，由最小二乘法得到重建信号并选择最大的s个元素所对应的原子,组成新的的支撑集F；
 (5)由初始余量，支撑集F，利用最小二乘法计算更新余量r；
 (6)如果满足迭代停止条件〖||r||〗_2 ≤ε(固定阈值)则退出循环，如果满足〖||r_k ||〗2 〖||r(k-1) ||〗_2则跟新步长倍数j = j + 1，跟新的步长为s = j × s转入步骤二继续迭代,否则更新支撑集，更新余量，转入步骤二继续迭代。
 为了说明SAMP算法的重建性能，利用MATLAB R2016a作为平台，分2次实验验证。第一次实验，假定信号是稀疏的并且稀疏度为10，定常信号256 × 1的列向量中随机选取10个随机数随机排列到定长信号中。固定阈值为0.1，观测矩阵为128×256高斯随机矩阵。初始信号和恢复信号之间的误差用二者相减的二范数表示。实验结果如(7)所示，误差为o(10 )。
 
第二次实验，测试对象为256×256lena图像。观测矩阵采用随机高斯矩阵，M和N表示观测矩阵的行数和列数，M/N表示压缩比，实验的压缩比为0.5，阈值为100，验的PSRN值为24.3。实验结果如图(8)
 
 
 
SAMP算法通过提供了严格的误差界限(固定阈值)且不需已知信号的稀疏度 K 就可以重构信号，固定阈值的选取同测量矩阵有着很大的联系。初始步长的选取目前任然是一个问题，SAMP算法的初始步长只需要小于稀疏度K，为了避免过度检测，如果不知道稀疏度,选择初始步长为1可以确保足够的精度。但是初始步长越小，算法运行的时间也就越大，这就造成了过度检测与运行时间的矛盾。经验表明，对于指数衰减的信号初始步长选的较小比较好，对于二进制稀疏信号选择步长较大比较好。如果算法能随着算法逐渐接近真实的稀疏度K值而逐渐缩小这样会提高算法的精度。另外，有学者提出了稀疏度自适应子空间追踪算法，该算法能较为精确完成信号重构，不需要设定步长，能够估计稀疏度大小，但是输入参数对算法影响较大。也有学者提出变步长的自适应步长的稀疏度估计的方法[12]，根据每次迭代得到的残差值，通过函数Ln=[L_(n-1)log⁡(γ)/log⁡((γ+0.9ε)) /2.1],其中γ=1-|(|r_2 |)|/||y_2 ||，来确定步长的大小，当步长远离稀疏度K时步长较大，当接近K时稀疏度较小。
 3.5基于小波变换的分块压缩感知
 根据压缩感知理论，图像重构时如果图像尺寸太大，为了维持一定的精度，测量矩阵所需要的观测值随之增加，造成了观测矩阵对于有限的存储空间显得十分巨大，而硬件难于满足需求。Lu Gan [13]提出了分块压缩感知12，该方法指出: 可以将原始图像分成一些大小相等的图像块，采用相同的观测矩阵单独对每个图像块进行观测，大大简化了计算复杂度，这种方法能解决大尺度图像实时传输的问题。分块图像的观测矩阵远远小于未分块的图像，降低了对硬件的需求，利用PC端将采集的样本重构恢复，并实现图像融合。
 本实验利用MATLAB R2016a作为平台，对每一个子块采用小波变换，保留每块的低频小波稀疏，对高频小波稀疏进行采样得到测量向量。13重构时利用正交匹配追踪(OMP) 算法对高频系数进行恢复，再进行小波反变换重构图像。
 实验步骤：
 选取大小为256 × 256的二维图像，将图像分为16块大小为64 × 64的图像；
 采用小波变换对图像进行稀疏化表示，利用50 × 64的观测矩阵进行观测得到观测向量；
 采用OMP算法对观测向量进行重构恢复；
 对重构的恢复矩阵进行小波反变换，得到重构图像。
 实验结果如图(9)所示
 
 
 
为了比较几种算法的恢复精度，在固定的稀疏度12下对上述算法进行实验，实验信号为定常信号256 × 1的列向量，观测随着观测矩阵行数M的数目变化对一维信号恢复正确率的变化，实验中对每次M的取值实验100次，计算恢复率，信号之间的误差达到10 即视为正确，M的取值为从稀疏度K每次增加十个直到观测矩阵列数256。实验结果如下图所示(10)所示
 
图(10)为四种算法随着测量值M取值的不同，恢复值百分数变化
 图中可以看出ROMP算法到达较高的精确度所需要的观测矩阵行数较其余的算法高，压缩量较低。SAMP算法在稀疏度未知的情况下，由较小的测量值得到恢复精确度较高。
 4.总结
 压缩感知利用信号稀疏的特性将原来基于奈奎斯特采样定理的信号采样过程转化为基于优化计算恢复信号的观测过程。有效缓解了高速采样实现的压力,减少了处理、存储和传输的成本,使得用低成本的传感器将模拟信息转化为数字信息成为可能，同时压缩感知使得信号的恢复率理论上比传统压缩算法更加精准。
 研究还存在如下问题：
 (1)对于一个稳定的优化算法，是否存在最优的观测矩阵使得观测值在达到一定的精度范围而观测数目控制在较少的范围。
 (2)如何设计有效的软硬件来应用压缩感知理论解决大量的实际问题,这方面的研究还不成熟。
 (3)含噪信号或采样过程中引入噪声时的信号重构问题也是难点所在。
 5.参考文献
 [1] 石光明, 刘丹华, 高大化等. 压缩感知理论及其研究进展[D]. 2009.
 [2] 李卓凡, 闫敬文. 压缩感知及应用[J]. 微计算机应用, 2010, 31(03): 12–16.
 [3] 李珅, 马彩文, 李艳等. 压缩感知重构算法综述[J]. 红外与激光工程, 2013, 42(S1): 225–232.
 [4] 吴赟, 丛琳. 压缩感知测量矩阵的研究[D]. 西安电子科技大学, 2012.
 [5] DONOHO D L. Compressed sensing[J]. IEEE Transactions on Information Theory, 2006, 52(4): 1289–1306.
 [6] NEEDELL D, VERSHYNIN R. Signal Recovery From Incomplete and Inaccurate Measurements Via Regularized Orthogonal Matching Pursuit[J]. Ieee Journal of Selected Topics in Signal Processing, 2010, 4(2): 310–316.
 [7] 杨真真, 杨震, 孙林慧. 信号压缩重构的正交匹配追踪类算法综述[J]. 信号处理, 2013, 29(04): 486–496.
 [8] TROPP J A, GILBERT A C. Signal Recovery From Random Measurements Via Orthogonal Matching Pursuit[J]. IEEE Transactions on Information Theory, 2007, 53(12): 4655–4666.
 [9] 汪浩然, 夏克文, 牛文佳. 分段正交匹配追踪(StOMP)算法改进研究[J]. 计算机工程与应用, 2017, 53(16): 55–61.
 [10] DO T T, GAN L, NGUYEN N等. Sparsity adaptive matching pursuit algorithm for practical compressed sensing[C]//2008 42nd Asilomar Conference on Signals, Systems and Computers. 2008: 581–587.
 [11] 杜秀丽, 胡兴, 顾斌斌等. 基于变步长的正则回溯SAMP压缩感知重构算法[J]. 计算机应用研究, 2018, 35(04): 1084–1087.
 [12] FU Y, LIU S, REN C. Adaptive Step-Size Matching Pursuit Algorithm for Practical Sparse Reconstruction[J]. Circuits Systems and Signal Processing, 2017, 36(6): 2275–2291.
 [13] 曹玉强, 柏森, 曹明武. 图像自适应分块的压缩感知采样算法[J]. 中国图象图形学报, 2016, 21(04): 416–424.

压缩感知常见稀疏基matlab实现代码
 
离散余弦变换矩阵
 
matlab有直接用的函数dctmtx。
 dctmtx(N)即何生成N×N的离散余弦正交变换矩阵。
 
离散傅里叶变换矩阵
 
matlab中同样有直接用的函数dftmtx。
 dftmtx(N)/sqrt(N)即可生成离散傅里叶变换矩阵。
 
正交小波变换矩阵
 
function ww=DWT(N)

[h,g]= wfilters('sym8','d');       %  分解低通和高通滤波器  

% N=256;                           %  矩阵维数(大小为2的整数幂次)
L=length(h);                       %  滤波器长度
rank_max=log2(N);                  %  最大层数
rank_min=double(int8(log2(L)))+1;  %  最小层数
ww=1;   %  预处理矩阵

%  矩阵构造
for jj=rank_min:rank_max
    
    nn=2^jj;
    
    %  构造向量
    p1_0=sparse([h,zeros(1,nn-L)]);
    p2_0=sparse([g,zeros(1,nn-L)]);
    
    %  向量圆周移位
    for ii=1:nn/2
        p1(ii,:)=circshift(p1_0',2*(ii-1))';
        p2(ii,:)=circshift(p2_0',2*(ii-1))';
    end
    
    %  构造正交矩阵
    w1=[p1;p2];
    mm=2^rank_max-length(w1);
    w=sparse([w1,zeros(length(w1),mm);zeros(mm,length(w1)),eye(mm,mm)]);
    ww=ww*w;
    
    clear p1;clear p2;
end

 
程序作者为沙威，香港大学电气电子工程学系，wsha@eee.hku.hk

 原文链接：http://blog.csdn.net/chulefei/article/details/51251916
 

 
 
 一、The Shannon-Nyquist Sampling Theorem
 
 问题：原始数据是连续函数，是否能用有限个采样百分之百重现原始数据？
 
 香农回答了这个问题：如果原始数据中最大频率为f，如果采样频率为2f，即每隔1/(2f)秒取一次样，则可完全恢复原始数据。
 
 
 
 
 陆吾生教授2010年的视频中给出了非常直观的解释：
 
 
 
 
 图a是采样和恢复过程，图b表示原始函数的傅里叶变换得到的频域分布，图c表示取样后的频域分布，可以看到是原始频域分布复制粘贴，且按采样频率位移，因为频率分布不能重叠，否则信息就会丢失，这个用公式表示就是下图的关系，所以可以很容易得到采样频率和原始频率的关系，图d是低通滤波器，图e是低通部分，即原始数据的频率分布。
 
 
 
 
 
 最终的恢复公式就是，可以看成是采样后的信号与sinc函数卷积的结果：
 
 
 
 
 其中sinc(x)部分是，它的图示曲线是：
 
 
 
 
 
 
 
 
 二、Sparse and Compressible Signals
 
 （我的理解）如果信号是由某种分布或者几种分布构成，则信号是可以被压缩的。
 
 比如离散余弦变换DCT和离散小波变换DWT，本质上即是预先设定的两组正交基，通过θ=C’x或θ=W‘x即可将原始向量解释为预设空间中的坐标。将θ中接近0的部分设为0，得到θ‘，则转换后的数据x' = Cθ'或x' = Wθ'，其中C'和W’分别为C和W的转置。
 
 
 
 
 
 
  
 
   [plain] 
   view plain
    copy 
   
 
   
 
    
   
   
  
 
 
 
 
 %Below is a MALAB code to generate matrix DCT.   
 function C = gen_dct(n)   
 alp = [sqrt(1/n) sqrt(2/n)*ones(1,n-1)];   
 ind = (1:2:(2*n-1))*pi/(2*n);   
 C = zeros(n,n);   
 for k = 1:n,   
 C(k,:) = alp(k)*cos((k-1)*ind);   
 end   
 
 
 
 
 
  
 
   [plain] 
   view plain
    copy 
   
 
   
 
    
   
   
  
 
 
 
 
 %Our second example is about an orthonormal basis based on discrete wavelet transform <span style="font-family: Arial, Helvetica, sans-serif; font-size: 12px;">DWT.</span>  
 % Suppose signal length n= 2^p  
 %for some integer p, then the basis matrix Wof size nby n  
 %associated with Daubechies’ orthogonal discrete wavelets can be readily constructed in MATLAB   
 %using Toolbox Uvi_Wave toolbox as follows:   
 % n -- signal length, must be a power of 2.   
 % L -- length of Daubechies wavelet, must be an even integer.   
 % W -- orthonormal DWT matrix of size n x n.   
 function W = gen_wave(n,L)   
 p = log2(n);   
 [h0,h1,f0,f1] = daub(L);   
 I = eye(n,n);   
 W = I;   
 for i = 1:n,   
 Ii = I(:,i);   
 W(:,i) = wt(Ii,h0,h1,p);   
 end   
 
 

 
 
 
 
 
 上面两图为余弦变换和小波变换的结果，50,100为对应的阈值，可以看出不同的正交基得到的结果不尽相同。
 
 陆吾生教授的视频里小波变换讲得非常好，想学习的一定要去看看。
 
 三、Sensing a Sparse Signal
 
 小波变换W，θ= W‘x，其实可以理解为用W中基底解释x的过程，x的稀疏性取决于W基底的”完备性“。这个可以理解为W是一个字典，W中的基底就是字典中的词，词越多你对原句的解释就可以越简洁，比如”上海“可以解释为”上“+”海“，也可以解释为”沪“，因此增加一个矩阵得到[I W]，可以得到更好的稀疏性。字典可以用其他已知的基底构成，比如单位矩阵I，离散余弦变换DCT，离散小波变换DWT，离散傅里叶变换DFT，Hadamard-Walsh矩阵等等。
 
 
 现在问题变为Dθ= x，D是n*l的字典，θ是l*1的解，x是输入向量，如何找到最稀疏的θ？
 
 Dθ = x其实是一个方程数小于变量数的线性方程组，所以这个方程的解一定是特解加通解的形式，特解是类似于广义逆的形式，θ1 = D‘inv(DD')x，通解需要对D进行SVD分解D=UΛV，U是n*n的特征向量矩阵，Λ是n*l特征值矩阵，V是l*l的特征向量矩阵，θ2 = Vn*δ，Vn是V后(l-n)列，即l*(l-n)的矩阵，δ是(l-n)*1的自由变量（随便是啥），问题就变成了如何选择δ确定的解空间中寻找最稀疏的θ1+θ2。
 
 最优化问题形式化为：
 
 
 
 
 其中0-norm表示θ中非0元素的个数，可是这个问题是NP-hard问题，不过可以转化为：
 
 
 
 
 1-norm就是绝对值的和，上面问题就可以转化为另一个不带绝对值的线性规划问题，这个问题可以用sedumi这个工具求解（具体过程请看陆教授的视频），得到的解比直接用DWT要好很多。
 
 四、Compressed Sensing
 
 
 
 按照香农的理论，假设采样频率为f=n，则一秒钟需要采样n个数据，从线性代数的角度来看，是将无限维原始数据映射到n维正交空间。DCT和DWT或者他们构成的字典提供了巧妙的正交基，使得新向量变得稀疏。现在有个想法，为什么不直接采样到好空间？换句话说，能否在采样时不是按照原始数据的频率，而是按照原始数据信息的频率。
 
 假设信号在Ψ空间中是r-sparse，转化为压缩信号θ：
 
 
 
 
 现在有另外一个空间Φ，采样x：
 
 
 
 
 Φ^{\hat}是Φ随机抽取m行的结果，y是m*1，Φ^{\hat}是m*n，x是n*1，m<n，所以这里是压缩采样。
 
 已知y重建x的过程就是最优化一个最稀疏的θ，然后用上面的公式求出x
 
 
 
 
 当m满足下面公式时，可以保证信号可以恢复：
 
 
 
 
 可以简化为：
 
 
 
 
 
 
 
 压缩感知的Matlab tool可以用CalTech的l1-MAGIC。
 
 Ψ的选择很重要，直接确定原始数据中哪些信息是noise哪些是真实信息，也就是决定信号是否可以压缩以及压缩后是否可以恢复，Φ可以是随机的矩阵。