CS-BP压缩感知的图像重建算法matlab仿真

Posted 2023-03-02 fpga和matlab

up目录

一、理论基础

二、核心程序

三、测试结果

---

一、理论基础

       压缩感知（Compressed Sensing，CS）指出只要信号是可压缩的或在某个变换域是稀疏的,那么就可以用一个与变换基不相关的观测矩阵将变换所得高维信号投影到一个低维空间上,然后通过求解一个优化问题就可以从这些少量的投影中以高概率重构出原信号。它包含两个特性，即不相关性和欠定性，压缩感知的压缩和重构正是靠这两个性质实现的。

在这一部分，需要明确几个名词的概念（顺序同符号说明表中符号顺序）：
原始信号：顾名思义，就是需要被处理的信号，可以是一维或二维信号；

稀疏信号：即原始信号投影后确定稀疏的信号。如何判断一个信号是否是稀疏信号，最简单的方法是验证该信号在某个变换域中是否只有极少个非零元素。比如，若长度为 N 的一维信号 x 中只有 K(K<<N) 个非零的元素，则称信号 x 是 K-稀疏 的；

稀疏基：又称正交基、稀疏矩阵，其主要作用是确保原始信号投影到某个域中的信号是稀疏的；

      观测矩阵：又称测量矩阵，控制采样得到的目标信号在采样过程中即被压缩，同时保证目标信号所含有效信息不丢失，能够由压缩采样值还原出目标信号；

       感知矩阵：又称传感矩阵、恢复矩阵，该矩阵其实是 稀疏基 与 观测矩阵 的乘积；

       观测值：又称观测信号，它是通过观测矩阵获得的值，其包含了信号中的有效信息，且数据量比原信号数据量小。

        从采样的角度来看，压缩感知和基于奈奎斯特采样定理的传统信号采集是两种不同形式的信号采集方式。（压缩感知打破了传统信号处理中对于奈奎斯特采样要求的限制）

       采样率：在压缩感知理论下，信号的采样率不再取决于信号的带宽，而是取决于信息在信号中的结构与内容（稀疏性）。关于采样率的计算方式，压缩感知是从少量离散测量数据恢复离散数字信号，其计算方式为采样率=测量值的大小/恢复信号的大小；而传统信号采集是从离散采样数据中恢复模拟信号（时序信号），采样率指的是一个采集频率，在我看来，这两者定义的采样率不具有可比性。（但从绝对值来看，压缩感知的采集数据量应该是小于或远小于传统采集）

       信号采集方式：传统采样理论是通过均匀采样（极少情况下也采用非均匀采样）获取数据；压缩感知则通过计算信号与一个观测函数之间的内积来获得观测数据（AX=b）；

       恢复信号形式：传统采样定理关注的对象是无限长的连续信号；压缩感知是有限维观测向量空间的向量（离散信号）；

        恢复信号方式：传统采样恢复是在Nyquist采样定理的基础上，通过采样数据的sinc函数线性内插获得（在不均匀采样下则是非线性的插值恢复），而压缩感知采用的是利用信号的稀疏性，从线性观测数据中通过求解一个非线性的优化问题来恢复信号的方法。

        压缩感知的核心思想：压缩和采样合并进行，并且测量值远小于传统采样方法的数据量，突破香农采样定理的瓶颈，使高分辨率的信号采集成为可能。

压缩感知理论主要包括三部分：

（1）信号的稀疏表示；

（2）设计测量矩阵，要在降低维数的同时保证原始信号x的信息损失最小；

（3）设计信号恢复算法，利用M个观测值无失真地恢复出长度为N的原始信号。

理论依据：

（1）设长度为N的信号X在某个正交基Ψ上是K-稀疏的（即含有k个非零值）；

（2）如果能找到一个与Ψ不相关（不相干）的观测基Φ；

（3）用观测基Φ观测原信号得到长度M的一维测量值M个观测值Y，K<M<<N；

（4）那么就可以利用最优化方法从观测值Y中高概率恢复X。

        当矩阵Φ满足RIP准则时。压缩感知理论能够通过对上式的逆问题先求解稀疏系数s，然后将稀疏度为K的信号x从M维的测量投影值y中正确地恢复出来。解码的最直接方法是通过l0范数（0-范数，也就是向量yˆ中非零元素的个数）下求解的最优化问题：

 

       从而得到稀疏系数s的估计s’。则原信号x’ = Ψs’。由于上式的求解是个NP难问题（在多项式时间内难以求解，甚至无法验证解的可靠性）。L1最小范数下在一定条件下和L0最小范数具有等价性，可得到相同的解。那么上式转化为L1最小范数下的最优化问题： 

  

       L 1范数最小化是通过用L1范数来近似0范数，取1而不取1/2,2/3或者其他值，是因为1范数最小化是凸优化问题，可以将求解过程转化成有一个线性规划问题。L1最小范数下最优化问题又称为基追踪(BP)，其常用实现算法有：内点法和梯度投影法。内点法速度慢，但得到的结果十分准确：而梯度投影法速度快，但没有内点法得到的结果准确 。 

二、核心程序

clc;
clear;
close all;
warning off;
addpath(genpath(pwd));


%------------ read in the image --------------
img=imread('lena.bmp');     % testing image
img=double(img);
[height,width]=size(img);

%------------ form the measurement matrix and base matrix ---------------
Phi=rand(floor(height/1.1),width);  % only keep one third of the original data  
Phi = Phi./repmat(sqrt(sum(Phi.^2,1)),[floor(height/1.1),1]); % normalize each column


mat_dct_1d=zeros(256,256);  % building the DCT basis (corresponding to each column)
for k=0:1:255 
    dct_1d=cos([0:1:255]'*k*pi/256);
    if k>0
        dct_1d=dct_1d-mean(dct_1d); 
    end;
    mat_dct_1d(:,k+1)=dct_1d/norm(dct_1d);
end


%--------- projection ---------
img_cs_1d=Phi*img;          % treat each column as a independent signal


%-------- recover using omp ------------
sparse_rec_1d=zeros(height,width);  

Theta_1d=Phi*mat_dct_1d;
    
for i=1:width
    column_rec=cs_bp(img_cs_1d(:,i),Theta_1d,height);
    sparse_rec_1d(:,i)=column_rec';           % sparse representation
end
img_rec_1d=mat_dct_1d*sparse_rec_1d;          % inverse transform


%------------ show the results --------------------
figure(1)
subplot(121)
imshow(img,[])
title('原图像');
subplot(122)
imshow(img_rec_1d,[])
title('压缩感知重建图像');

up93

三、测试结果