利用MATLAB实现白话压缩感知

pulbieup

1 V/ d+ W3 \4 R. q7 D: }8 _! z8 \压缩感知介绍
压缩感知（Compressive Sensing，CS），有时也叫成Compressive Sampling。相对于传统的奈奎斯特采样定理——要求采样频率必须是信号最高频率的两倍或两倍以上（这就要求信号是带限信号，通常在采样前使用低通滤波器使信号带限），压缩感知则利用数据的冗余特性，只采集少量的样本还原原始数据。
  A6 Q7 X7 n5 \
这所谓的冗余特性，借助MLSS2014马毅老师的课件上的例子来说明，

* x+ T# V8 C9 H3 H! V
  T, h: t/ y' d: ^# j$ z
) N# ]8 A6 c3 h/ k( M5 C, H因为自然界的数据都存在局部低维结构、周期性、对称性等，因此，传统的固定采样率的采样方法必然存在信息冗余。由于信息冗余的存在，我们就知道有数据压缩那么一门学科。既然这样，为什么要把冗余的数据也采进来，再进行压缩呢，能不能不把冗余的数据采集进来？

压缩感知的思路就是：在采集的过程中就对数据进行了压缩，而且这种压缩能保证不失真（低失真）的恢复原始数据，这与传统的先2倍频率采集信号→存储→再压缩的方式不同，可以降低采集信号的存储空间和计算量。
6 V  b1 x6 |  W8 O! U' [! q6 K+ F
好了，那么就来了解一下压缩感知的具体模型。

" s) h! Q- _8 k  T! s
" [/ T0 W* q7 E6 f# @/ Z1. 稀疏表示

使用压缩感知理论首先要求信号能表示为稀疏信号，如x=[1 0 0 0 1 0]，其中只有2个1，可认为是稀疏的。我们将信号通过一个矩阵映射到稀疏空间，
# g  H2 i, r, M. }* T. m
设信号x为N维，即，则为NxN维稀疏表达矩阵，s即是将x进行稀疏表示后的Nx1维向量，其中大部分元素值为0。稀疏表示的原理就是通过线性空间映射，将信号在稀疏空间进行表示。

比如，信号

在时域是非稀疏的，但做傅里叶变换表示成频域后，只有少数几个点为非0（如下图）。则该信号的时域空间为非稀疏，频域空间为稀疏空间，组成的矩阵。一般为正交矩阵。
+ N( e  ]# v$ ?" J
9 h. ^+ j/ M( d! u0 t/ ]6 Q

: w5 ^2 h0 ]' p4 [: ~若稀疏表示后的结果s中只有k个值不为0，则称x的稀疏表示为k-Sparse。上图对x的频域稀疏表示就是2-Sparse。
5 Y" F" M, n. M- I# V0 K( _8 k
4 m  |5 d3 i. x, v+ r; G
2. 感知测量
0 M( f2 ~! l- k) H( \
1 n, a5 |* y0 N
7 o4 M# Y9 Q% `1 l压缩感知的目的是在采集信号时就对数据进行压缩，大牛们的思路集中到了数据采集上——既然要压缩，还不如就从大量的传感器中只使用其中很少的一部分传感器，采集少量的冗余度低的数据。这就是感知测量的通俗的说法，用表达式表示
  o/ d. h: [! K, n+ n; a6 S2 F
5 U! x/ i* M! z其中的x就是稀疏表示中的信号，为MxN维的感知矩阵（M表示测量信号的维度），y则表示M（M远小于N才有意义）个传感器的直接测量，因此维度为Mx1。
: n; |: @# T0 }
; S# O/ d2 |9 z. t' A/ i6 H& C% |  j将稀疏表示过程和感知测量过程综合起来：
- L, r% }2 c$ O7 ]  @, o
$ l, V* c7 `- p+ X
5 g; f/ O. C2 n3 i2 k6 g7 b
% v# z+ T! K& p
数学描述

, I7 r; T- X. }: `: N" c
对于压缩感知模型，其中每个量的维度一定要了解（通过维度的变化来理解压缩感知很有效）：
# ^+ L! p& d* [( B. V9 ]


$ e4 M& j) z; N: @- m2 _* A从感知测量中知道：M就是测量的维度（M远小于N）。
1 }& w# |. L. v5 P0 [' Z2 I# c& S
2 s4 K9 T. n8 e6 t. m压缩感知的原信号恢复问题描述为：
) E9 t1 c& D$ h- P. S6 C* q
由已知条件：
0 n- H" D  _) d& T) W
3 }' r0 T9 l8 l, U(1) 测量值y，且，其中e为噪声引入
/ i- x. G5 B; n& |
  ^1 y8 P; k2 m2 I5 A(2) s为k-Sparse信号（k未知）
) e& w. X  v# J6 B( w7 W" w1 D
求解目标：k尽可能小的稀疏表示信号s及对应的。

用数学形式描述为：
/ d9 V. t9 A0 h3 ^) A
! @) ]2 F5 D$ m. w" Ue可以看成告诉随机噪声，e~N(0,δ2)。

# K6 o9 m: w" r- j即是要求s使s的0范数（非0值的个数）最小，但0范数优化问题是很难求解的，于是一帮大牛就证明求解1范数也能逼近和上面相同的效果，而求解2范数及其更高的范数则结果相差越来越大（有些人在研究介于0范数与1范数之间的范数求解方法）。因此可转化为1范数求解：
7 D0 u. I- s/ P7 K
由拉格朗日乘子法，上面的最优问题可转化成：

' S4 l2 m) M+ y上面的最小值求解问题就可以直接通过凸优化解得结果了。

0 e+ @/ f: n" j* D+ \2 U0 A
  h5 f) i2 v4 m0 N5 b1 r程序分析
% F4 y& y+ h$ ?2 s" ?, x

先下载CVX或spams工具箱，下面的matlab程序分别使用了时域稀疏信号和频域稀疏信号进行测试，两种不同在于时域稀疏信号不用稀疏表达矩阵（因此稀疏表达矩阵使用单位阵即可），而频域稀疏信号则需要先通过稀疏表达矩阵将信号在频域进行稀疏表示，再压缩感知后进行恢复时也要进行反FFT变换到时域。
) I. v) B8 h( P
关于M的选取：测量数M>=K*log(N/K),K是稀疏度,N信号长度,可以近乎完全重构
& V; [0 z% n. m( N5 T& w* B- f/ ?
8 Z- h4 L2 e& w* f% _clc
clear all
* e1 F& X% }& `# n; Eclose all

+ b) h& d  L8 f& l4 `%% 产生原始信号及相关参数
n = 512;                          % 信号长度
s = 25;                           % 稀疏度（从下面知道，分时域和频域的情况）
m = 5*s;                          % 测量长度 M>=S*log(N/S)
! ]. J1 ~: v! Y7 zfreq_sparse = 0;                  % 信号频域稀疏则为1

: I; M$ w! k0 m, j- g. ]if freq_sparse
    t = [0: n-1]';
    f = cos(2*pi/256*t) + sin(2*pi/128*t);   % 产生频域稀疏的信号
else
& u# H. ^/ ?( G6 d  C" q! K, ?( a4 w    tmp = randperm(n);   
) q0 S2 w/ _' ]& u1 j    f = zeros(n,1);
/ [* v+ D  `6 F: p    f(tmp(1:s)) = randn(s,1);     % 产生时域稀疏的信号
) {! `2 \* M* O; a+ y$ r  k& U- ^end

( x. D) W! V& o% `" {3 X%% 产生感知矩阵和稀疏表示矩阵
$ ~* s% G$ f" ]$ X% F' _Phi = sqrt(1/m) * randn(m,n);     % 感知矩阵（测量矩阵）
% A = get_A_fourier(n, m);

) J2 ^0 d  x$ hy = Phi * f;                      % 通过感知矩阵获得测量值
" o# v$ O$ F" H% r& Q+ H. G                                  % Psi 将信号变换到稀疏域
4 r1 q9 l  |2 y7 I: E" O9 C+ \& Z2 `if freq_sparse                    % 信号频域可以稀疏表示
    Psi = inv(fft(eye(n,n)));     % 傅里叶正变换，频域稀疏正交基（稀疏表示矩阵）
else                              % 信号时域可以稀疏表示
    Psi = eye(n, n);              % 时域稀疏正交基
end
5 v+ }* I, Z4 JA = Phi * Psi;                    % 恢复矩阵 A = Phi * Psi
3 g/ ^! f. s) s' }
3 K2 y/ S# N3 C%% 优化方法1：使用CVX工具进行凸优化
addpath('../../cvx-w64/cvx/');
cvx_startup;                            % 设置cvx环境
cvx_begin
) _, U" j5 J0 t% n% O" V5 t  s    variable xp(n) complex;               % 如果xp是复数，则添加complex，否则不加
    minimize (norm(xp, 1));
6 [7 i" k7 R  B    subject to
        A * xp == y;      
cvx_end
! n$ D8 d* m; r2 m) P
%% 优化方法2：使用spams工具箱进行优化
% addpath('spams-matlab/build');
% param.L = 100;
% param.eps = 0.001;
8 w! H5 ]4 ?# W7 ?. ]% param.numThreads = -1;
%
; I: p/ [+ A0 Q4 A% 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
; h8 @) g( J4 j- `
end
plot(f+noise);
hold on
( _& m8 h2 m+ @* u1 }. \: _! f/ [plot(real(xp), 'r.');
legend('Original', 'Recovered');
+ Q+ p3 }8 p! f* k" M: K

• 设程序中的freq_sparse = 0时，观察时域稀疏信号的恢复结果为
  

  G  J2 R7 L! m8 [, M3 |) S
在恢复时域稀疏信号时，所使用的测量矩阵Phi为初始化的随机阵，因为本身时域就稀疏，而算法直接在时域进行恢复，所以稀疏表达矩阵就是单位阵Psi=E。上面的时域稀疏恢复结果显得没有误差是因为没有给原始信号添加噪声。
8 P! r* F" d, F* e5 y' L$ m
7 c; a4 r2 O" ]2. 设程序中的freq_sparse = 1时，观察频域稀疏信号的恢复结果为
( w1 C$ ~( d9 u
$ n, C7 U0 \  ?0 x

% c; G( h0 [- `5 V, i在恢复时域稀疏信号时，所使用的测量矩阵Phi为初始化的随机阵，因为信号只在频域稀疏，所以稀疏变换矩阵为傅里叶变换基，所以稀疏表达矩阵就是Psi = inv(fft(eye(n,n)))。同理，上面的频域稀疏恢复结果显得没有误差是因为没有给原始信号添加噪声。
+ p/ v4 M* L6 w! Q, v: p
7 A1 X6 P7 J- U/ P' F( T+ ^- l3. 上面都是没有添加噪声的算法结果，我们给信号添加一些噪声，以时域信号为例，

0 \0 v" H0 T) G4 q, M% h' j) F2 dif freq_sparse
3 z! W( w5 C; D1 q* S+ h$ p    t = [0: n-1]';
    f = cos(2*pi/256*t) + sin(2*pi/128*t);   % 产生频域稀疏的信号
" s/ ]1 D  o, \2 t  |# Oelse
/ x4 C2 |3 N& U+ V( |9 u    tmp = randperm(n);   
    f = zeros(n,1);
" P. K9 {, g% h    f(tmp(1:s)) = randn(s,1);     % 产生时域稀疏的信号
end

6 r/ H/ C  R, Z/ ?) }# tnoise = random('norm', 0, 0.01, [n 1]);
% z3 b" P8 G/ H1 f5 m; ~2 Xf = f + noise;                    % 添加噪声
; m) S9 ?+ m* J+ ]
%% Remaining code not changed


从下图的恢复结果看，已经能看到一些恢复误差了，
! F/ r- B7 i2 v  C  M, [; O) \
/ I0 }- J" j+ k4 P/ \

ExxNEN

利用MATLAB实现白话压缩感知