壓縮感知介紹

壓縮感知（Compressive Sensing，CS），有時也叫成Compressive Sampling。相對於傳統的奈奎斯特取樣定理——要求取樣頻率必須是訊號最高頻率的兩倍或兩倍以上（這就要求訊號是帶限訊號，通常在取樣前使用低通濾波器使訊號帶限），壓縮感知則利用資料的冗餘特性，只採集少量的樣本還原原始資料。

這所謂的冗餘特性，藉助MLSS2014馬毅老師的課件上的例子來說明，

因為自然界的資料都存在區域性低維結構、週期性、對稱性等，因此，傳統的固定取樣率的取樣方法必然存在資訊冗餘。由於資訊冗餘的存在，我們就知道有資料壓縮那麼一門學科。既然這樣，為什麼要把冗餘的資料也採進來，再進行壓縮呢，能不能不把冗餘的資料採集進來？

壓縮感知的思路就是：在採集的過程中就對資料進行了壓縮，而且這種壓縮能保證不失真（低失真）的恢復原始資料，這與傳統的先2倍頻率採集訊號→儲存→再壓縮的方式不同，可以降低採集訊號的儲存空間和計算量。

好了，那麼就來了解一下壓縮感知的具體模型。

1. 稀疏表示

使用壓縮感知理論首先要求訊號能表示為稀疏訊號，如x=[1 0 0 0 1 0]，其中只有2個1，可認為是稀疏的。我們將訊號通過一個矩陣對映到稀疏空間，

設訊號x為N維，即，則為NxN維稀疏表達矩陣，s即是將x進行稀疏表示後的Nx1維向量，其中大部分元素值為0。稀疏表示的原理就是通過線性空間對映，將訊號在稀疏空間進行表示。

比如，訊號

在時域是非稀疏的，但做傅立葉變換表示成頻域後，只有少數幾個點為非0（如下圖）。則該訊號的時域空間為非稀疏，頻域空間為稀疏空間，

若稀疏表示後的結果s中只有k個值不為0，則稱x的稀疏表示為k-Sparse。上圖對x的頻域稀疏表示就是2-Sparse。

2. 感知測量

壓縮感知的目的是在採集訊號時就對資料進行壓縮，大牛們的思路集中到了資料採集上——既然要壓縮，還不如就從大量的感測器中只使用其中很少的一部分感測器，採集少量的冗餘度低的資料。這就是感知測量的通俗的說法，用表示式表示

其中的x就是稀疏表示中的訊號，為MxN維的感知矩陣（M表示測量訊號的維度），y則表示M（M遠小於N才有意義）個感測器的直接測量，因此維度為Mx1。

將稀疏表示過程和感知測量過程綜合起來：

數學描述

對於壓縮感知模型，其中每個量的維度一定要了解（通過維度的變化來理解壓縮感知很有效）：

從感知測量中知道：M就是測量的維度（M遠小於N）。

壓縮感知的原訊號恢復問題描述為：

由已知條件：

(1) 測量值y，且，其中e為噪聲引入

(2) s為k-Sparse訊號（k未知）

求解目標：k儘可能小的稀疏表示訊號s及對應的。

用數學形式描述為：

e可以看成告訴隨機噪聲，e~N(0,δ2)。

即是要求s使s的0範數（非0值的個數）最小，但0範數優化問題是很難求解的，於是一幫大牛就證明求解1範數也能逼近和上面相同的效果，而求解2範數及其更高的範數則結果相差越來越大（有些人在研究介於0範數與1範數之間的範數求解方法）。因此可轉化為1範數求解：

由拉格朗日乘子法，上面的最優問題可轉化成：

上面的最小值求解問題就可以直接通過凸優化解得結果了。

程式分析

先下載CVX或spams工具箱，下面的matlab程式分別使用了時域稀疏訊號和頻域稀疏訊號進行測試，兩種不同在於時域稀疏訊號不用稀疏表達矩陣（因此稀疏表達矩陣使用單位陣即可），而頻域稀疏訊號則需要先通過稀疏表達矩陣將訊號在頻域進行稀疏表示，再壓縮感知後進行恢復時也要進行反FFT變換到時域。

關於M的選取：測量數M>=K*log(N/K),K是稀疏度,N訊號長度,可以近乎完全重構

clc
clear all
close all

%% 產生原始訊號及相關引數
n = 512;                          % 訊號長度
s = 25;                           % 稀疏度（從下面知道，分時域和頻域的情況）
m = 5*s;                          % 測量長度 M>=S*log(N/S)
freq_sparse = 0;                  % 訊號頻域稀疏則為1

if freq_sparse
    t = [0: n-1]';
    f = cos(2*pi/256*t) + sin(2*pi/128*t);   % 產生頻域稀疏的訊號
else
    tmp = randperm(n);    
    f = zeros(n,1);
    f(tmp(1:s)) = randn(s,1);     % 產生時域稀疏的訊號
end

%% 產生感知矩陣和稀疏表示矩陣
Phi = sqrt(1/m) * randn(m,n);     % 感知矩陣（測量矩陣）
% A = get_A_fourier(n, m);

y = Phi * f;                      % 通過感知矩陣獲得測量值
                                  % Psi 將訊號變換到稀疏域
if freq_sparse                    % 訊號頻域可以稀疏表示
    Psi = inv(fft(eye(n,n)));     % 傅立葉正變換，頻域稀疏正交基（稀疏表示矩陣）
else                              % 訊號時域可以稀疏表示
    Psi = eye(n, n);              % 時域稀疏正交基
end
A = Phi * Psi;                    % 恢復矩陣 A = Phi * Psi

%% 優化方法1：使用CVX工具進行凸優化
addpath('../../cvx-w64/cvx/');
cvx_startup;                            % 設定cvx環境
cvx_begin
    variable xp(n) complex;               % 如果xp是複數，則新增complex，否則不加
    minimize (norm(xp, 1));
    subject to
        A * xp == y;      
cvx_end

%% 優化方法2：使用spams工具箱進行優化
% addpath('spams-matlab/build');
% param.L = 100;
% param.eps = 0.001;
% param.numThreads = -1;
% 
% A=A./repmat(sqrt(sum(A.^2)),[size(A,1) 1]);
% xp = mexOMP(y, A, param);       % 正交匹配追蹤法Orthogonal Matching Pursuit

%% 對比原訊號和
if freq_sparse
    xp = real(ifft(full(xp)));           % 傅立葉逆變換
else

end
plot(f+noise);
hold on
plot(real(xp), 'r.');
legend('Original', 'Recovered');

1. 設程式中的freq_sparse = 0時，觀察時域稀疏訊號的恢復結果為

   

   在恢復時域稀疏訊號時，所使用的測量矩陣Phi為初始化的隨機陣，因為本身時域就稀疏，而演算法直接在時域進行恢復，所以稀疏表達矩陣就是單位陣Psi=E。上面的時域稀疏恢復結果顯得沒有誤差是因為沒有給原始訊號新增噪聲。

2. 設程式中的freq_sparse = 1時，觀察頻域稀疏訊號的恢復結果為

   

   在恢復時域稀疏訊號時，所使用的測量矩陣Phi為初始化的隨機陣，因為訊號只在頻域稀疏，所以稀疏變換矩陣為傅立葉變換基，所以稀疏表達矩陣就是Psi = inv(fft(eye(n,n)))。同理，上面的頻域稀疏恢復結果顯得沒有誤差是因為沒有給原始訊號新增噪聲。

3. 上面都是沒有新增噪聲的演算法結果，我們給訊號新增一些噪聲，以時域訊號為例，

   if freq_sparse
       t = [0: n-1]';
       f = cos(2*pi/256*t) + sin(2*pi/128*t);   % 產生頻域稀疏的訊號
   else
       tmp = randperm(n);    
       f = zeros(n,1);
       f(tmp(1:s)) = randn(s,1);     % 產生時域稀疏的訊號
   end
   
   noise = random('norm', 0, 0.01, [n 1]);
   f = f + noise;                    % 新增噪聲
   
   %% Remaining code not changed
   

   從下圖的恢復結果看，已經能看到一些恢復誤差了，