从稀疏表示到压缩感知（下）

原From <http://xiahouzuoxin.github.io/notes/html/%E7%99%BD%E8%AF%9D%E5%8E%8B%E7%BC%A9%E6%84%9F%E7%9F%A5.html>

python代码自己添加

  

1. 稀疏表示

使用压缩感知理论首先要求信号能表示为稀疏信号，如x=[1 0 0 0 1 0]，其中只有2个1，可认为是稀疏的。我们将信号通过一个矩阵映射到稀疏空间，

设信号x为N维，即

，则

为NxN维稀疏表达矩阵，s即是将x进行稀疏表示后的Nx1维向量，其中大部分元素值为0。稀疏表示的原理就是通过线性空间映射，将信号在稀疏空间进行表示。

 

  

1.1傅里叶变换python下的演示

原型：

import numpy as np

import matplotlib.pyplot as plt

  

# Compute the x and y coordinates for points on a sine curve

#x = np.arange(0, 3 * np.pi, 0.1)

#y = np.sin(x)

t = np.arange(-300*np.pi,300*np.pi,0.1)

x1=np.cos(2*np.pi*t/256)+np.sin(2*np.pi*t/128)

  

# Plot the points using matplotlib

plt.plot(t, x1)

plt.show() # You must call plt.show() to make graphics appear.

  

  

结果如图：

  

使用傅里叶对x1进行变换：

  

备注：

在numpy下傅里叶变换的函数

numpy.fft

for definition of the DFT and conventions used.

ifft

The inverse of fft.

fft2

The two-dimensional FFT.

fftn

The n-dimensional FFT.

rfftn

The n-dimensional FFT of real input.

fftfreq

Frequency bins for given FFT parameters.

  

From <http://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.fft.html#numpy.fft.fft>

 

  

#coding=utf-8

import numpy as np

import matplotlib.pyplot as plt

  

# Compute the x and y coordinates for points on a sine curve

#x = np.arange(0, 3 * np.pi, 0.1)

#y = np.sin(x)

t = np.arange(-300*np.pi,300*np.pi,0.1)

x1=np.cos(2*np.pi*t/256)+np.sin(2*np.pi*t/128)

  

#对x1进行傅里叶变换，反变换使用ifft

x2=np.fft.fft(x1)

# Plot the points using matplotlib

plt.plot(t, x2.real,t,x2.imag)

plt.show() # You must call plt.show() to make graphics appear.

结果：

可以看出来在时域是非稀疏的，但做傅里叶变换表示成频域后，只有少数几个点为非0（如下图）。则该信号的时域空间为非稀疏，频域空间为稀疏空间，组成的矩阵。一般为正交矩阵。

  

若稀疏表示后的结果s中只有k个值不为0，则称x的稀疏表示为k-Sparse。上图对x的频域稀疏表示就是2-Sparse。

———————————————————————————————————————————————————————————————————————

2. 感知测量

压缩感知的目的是在采集信号时就对数据进行压缩，大牛们的思路集中到了数据采集上——既然要压缩，还不如就从大量的传感器中只使用其中很少的一部分传感器，采集少量的冗余度低的数据。这就是感知测量的通俗的说法，用表达式表示

其中的x就是稀疏表示中的信号， 为MxN维的感知矩阵（M表示测量信号的维度），y则表示M（M远小于N才有意义）个传感器的直接测量，因此维度为Mx1。

将稀疏表示过程和感知测量过程综合起来：

  

———————————————————————————————————————————————————

数学描述

  

对于压缩感知模型，其中每个量的维度一定要了解（通过维度的变化来理解压缩感知很有效）：

yax

从感知测量中知道：M就是测量的维度（M远小于N）。

压缩感知的原信号恢复问题描述为：

由已知条件：

• 测量值y，且

  

  ，其中e为噪声引入

• 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);

hold on

plot(real(xp), 'r.');

legend('Original', 'Recovered');

  

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

  

  time
  在恢复时域稀疏信号时，所使用的测量矩阵Phi为初始化的随机阵，因为本身时域就稀疏，而算法直接在时域进行恢复，所以稀疏表达矩阵就是单位阵Psi=E。上面的时域稀疏恢复结果显得没有误差是因为没有给原始信号添加噪声。

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

  

  freq
  在恢复时域稀疏信号时，所使用的测量矩阵Phi为初始化的随机阵，因为信号只在频域稀疏，所以稀疏变换矩阵为傅里叶变换基，所以稀疏表达矩阵就是Psi = inv(fft(eye(n,n)))。同理，上面的频域稀疏恢复结果显得没有误差是因为没有给原始信号添加噪声。

• 上面都是没有添加噪声的算法结果，我们给信号添加一些噪声，以时域信号为例，

    

% 原信号增加噪声

noise = random('norm', 0, 0.01, [n 1]); % 噪声

f = f+ noise; % 信号添加噪声

  

% ....

% 以下是含噪声的优化，添加噪声项e

e = random('norm', 0, 0.01, [m 1]); % 噪声

addpath('../../cvx-w64/cvx/');

cvx_startup; % 设置cvx环境

cvx_begin

variable xp(n) complex; % 如果xp是复数，则添加complex，否则不加

minimize (norm(xp, 1));

subject to

A * xp + e == y; % 含噪声的优化

cvx_end

  

%% Remaining code not changed

  

从下图的恢复结果看，已经能看到一些恢复误差了