压缩感知和稀疏信号处理

一、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的转置。

%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

%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哪些是真实信息，也就是决定信号是否可以压缩以及压缩后是否可以恢复，Φ可以是随机的矩阵。