压缩感知的实现（含matlab代码）_压缩感知可以直接对矩阵做吗 ...

相关文章推荐

飘逸的黄瓜 · 心智圖在教育中的7種用法 - Edraw· 2 月前 ·

忧郁的烤土司 · 薇芯4470166 ...· 4 月前 ·

沉着的大白菜 · Curbside Recycling | ...· 6 月前 ·

腹黑的蚂蚁 · 【华工科技：预计800G光模块出货量将在 ...· 1 年前 ·

私奔的卤蛋 · 获取共享Outlook日历事件C#-腾讯云开 ...· 1 年前 ·

---------------------------------------------------------------------------------------------------------------------------------------------

（1）刘海英. 基于压缩感知理论的高光谱图像重建和超分辨成像技术研究[D]. 西安电子科技大学, 2014.

（2）压缩感知的常见稀疏基名称及离散傅里叶变换基

（3）压缩感知的常见测量矩阵

（1）压缩感知重构算法之正交匹配追踪(OMP)

----------------------------------------------------------------------------------------------------------------------------------------------

压缩感知（Compressive Sensing，CS）。相对于传统的奈奎斯特采样定理——要求采样频率必须是信号最高频率的两倍或两倍以上（这就要求信号是带限信号，通常在采样前使用低通滤波器使信号带限），压缩感知则利用数据的冗余特性，只采集少量的样本还原原始数据。

一句话总结我理解的压缩感知实现方法：以被重建信号在某个变换域上稀疏作为先验信息，用测量矩阵观测被测信号，由观测值结合重建算法重建出完整的被测信号。

在具体应用时，我们必须解决 CS 理论的三大关键问题：

目标信号的稀疏表示。寻找使得目标信号 f 变换到其上尽可能稀疏的变换域Ψ ，即信号稀疏表示问题；
测量矩阵的构建。测量矩阵是 CS 理论采样的实现部分。通过测量矩阵控制的采样使得目标信号 f在采样过程中即被压缩，同时保证目标信号所含有效信息不丢失，能够由压缩采样值还原出目标信号；
重建算法的设计。重建算法是从采样值求解最优化问题寻找到目标信号最优解。重建算法的准确性、高效性和稳定性是其设计的关键。

对于目标信号的稀疏表示问题，常见的稀疏基有离散余弦变换基（DCT）和快速傅立叶变换基（FFT）等。

对于测量矩阵，常见的有高斯随机矩阵、部分哈达玛矩阵等。

对于重建算法，常见的有L1范数、正交匹配追踪算法（OMP）等。

对于原理部分，相关文献、博客等资源相当多，本文不在这里赘述，详情可以参考本文开头引用内容。

本文分别以稀疏基有离散余弦变换基（DCT）和快速傅立叶变换基（FFT）做为稀疏基，高斯随机矩阵、部分哈达玛矩阵为测量矩阵，L1范数、正交匹配追踪算法（OMP）为重建算法进行压缩感知算法实现。

本文以f = cos(2*pi/256*t) + sin(2*pi/128*t)做为原信号，取原信号f的20%做为输入进行压缩感知重建。

注意：本文main.m中L1范数求解方法使用了CVX工具箱， ~~CVX工具箱安装方法参考（ CVX工具包（for matlab））~~ ，链接已失效，具体安装包自行百度。

main.m

%   该程序用于验证压缩感知理论（包含了L1最小范数求解和OMP求解）
clear all; close all;
%% 产生信号
choice_transform = 1;      % 选择正交基，1为选择DCT变换，0为选择FFT变换
choice_Phi = 0;         %选择测量矩阵，1为部分哈达玛矩阵，0为高斯随机矩阵
%-----------------------利用三角函数生成频域或DCT域离散信号--------------------------
n = 512;
t = [0: n-1];
f = cos(2*pi/256*t) + sin(2*pi/128*t);   % 产生频域稀疏的信号
%-------------------------------信号降采样率-----------------------
n = length(f);
a = 0.2;            %    取原信号的 a%
m = double(int32(a*n));
%--------------------------------------画信号图--------------------------------------
switch choice_transform
    case 1
        ft = dct(f);
        disp('ft = dct(f)')
    case 0
        ft = fft(f);
        disp('ft = fft(f)')
disp(['信号稀疏度：',num2str(length(find((abs(ft))>0.1)))])
figure('name', 'A Tone Time and Frequency Plot');
subplot(2, 1, 1);
plot(f);
xlabel('Time (s)'); 
% ylabel('f(t)');
subplot(2, 1, 2);
switch choice_transform
    case 1
        plot(ft)
        disp('plot(ft)')
    case 0
        plot(abs(ft));
        disp('plot(abs(ft))')
xlabel('Frequency (Hz)'); 
% ylabel('DCT(f(t))');
%% 产生感知矩阵和稀疏表示矩阵
%--------------------------利用感知矩阵生成测量值---------------------
switch choice_Phi
    case 1
        Phi = PartHadamardMtx(m,n);       % 感知矩阵（测量矩阵）    部分哈达玛矩阵
    case 0
        Phi = sqrt(1/m) * randn(m,n);     % 感知矩阵（测量矩阵）   高斯随机矩阵
% Phi =  randn(m,n);    %randn 生成标准正态分布的伪随机数（均值为0，方差为1）
% Phi = rand(m,n);    % rand 生成均匀分布的伪随机数。分布在（0~1）之间
f2 = (Phi * f')';                 % 通过感知矩阵获得测量值
% f2 = f(1:2:n);
switch choice_transform
    case 1
        Psi = dct(eye(n,n));            %离散余弦变换正交基 代码亦可写为Psi = dctmtx(n);
        disp('Psi = dct(eye(n,n));')
    case 0
        Psi = inv(fft(eye(n,n)));     % 傅里叶正变换，频域稀疏正交基（稀疏表示矩阵）
        disp('Psi = inv(fft(eye(n,n)));')
A = Phi * Psi;                    % 恢复矩阵 A = Phi * Psi
%%             重建信号
%---------------------使用CVX工具求解L1范数最小值-----------------
cvx_begin;
    variable x(n) complex;
%     variable x(n) ;
    minimize( norm(x,1) );
    subject to
      A*x == f2' ;
cvx_end;
figure;subplot(2,1,2)
switch choice_transform
    case 1
        plot(real(x));
        disp('plot(real(x))')
    case 0
        plot(abs(x));
        disp(' plot(abs(x))')
title('Using L1 Norm（Frequency Domain）');
%  ylabel('DCT(f(t))'); xlabel('Frequency (Hz)');
switch choice_transform
    case 1
        sig = dct(real(x));
        disp('sig = dct(real(x))')
    case 0
        sig = real(ifft(full(x)));
        disp(' sig = real(ifft(full(x)))')
subplot(2,1,1);
plot(f)
hold on;plot(sig);hold off
title('Using L1 Norm (Time Domain)');
% ylabel('f(t)'); xlabel('Time (s)');
legend('Original','Recovery')
%-----------------------------使用OMP算法重建--------------------------
for K = 1:100
    theta = CS_OMP(f2,A,K);
    %     figure;plot(dct(theta));title(['K=',num2str(K)])
    switch choice_transform
        case 1
            re(K) = norm(f'-(dct(theta)));
        case 0
            re(K) = norm(f'-real(ifft(full(theta))));
theta = CS_OMP(f2,A,find(re==min(min(re))));
disp(['最佳稀疏度K=',num2str(find(re==min(min(re))))]);
% theta = CS_OMP(f2,A,10);
figure;subplot(2,1,2);
switch choice_transform
    case 1
        plot(theta);
        disp('plot(theta)')
    case 0
        plot(abs(theta));
        disp('plot(abs(theta))')
title(['Using OMP(Frequence Domain)  K=',num2str(find(re==min(min(re))))])
switch choice_transform
    case 1
        sig2 = dct(theta);
        disp('sig2 = dct(theta)')
    case 0
        sig2 = real(ifft(full(theta)));
        disp('sig2 = real(ifft(full(theta)))')
subplot(2,1,1);plot(f);hold on;
plot(sig2)
hold off;
title(['Using OMP(Time Domain)  K=',num2str(find(re==min(min(re))))]);
legend('Original','Recovery')
部分哈达玛矩阵：PartHadamardMtx.m 
function [ Phi ] = PartHadamardMtx( M,N )  
%PartHadamardMtx Summary of this function goes here  
%   Generate part Hadamard matrix   
%   M -- RowNumber  
%   N -- ColumnNumber  
%   Phi -- The part Hadamard matrix  
% 来源http://blog.csdn.net/jbb0523/article/details/44700735 
%% parameter initialization  
%Because the MATLAB function hadamard handles only the cases where n, n/12,  
%or n/20 is a power of 2  
    L_t = max(M,N);%Maybe L_t does not meet requirement of function hadamard  
    L_t1 = (12 - mod(L_t,12)) + L_t;  
    L_t2 = (20 - mod(L_t,20)) + L_t;   
    L_t3 = 2^ceil(log2(L_t));  
    L = min([L_t1,L_t2,L_t3]);%Get the minimum L  
%% Generate part Hadamard matrix     
    Phi = [];  
    Phi_t = hadamard(L);  
    RowIndex = randperm(L);  
    Phi_t_r = Phi_t(RowIndex(1:M),:);  
    ColIndex = randperm(L);  
    Phi = Phi_t_r(:,ColIndex(1:N));  
正交匹配追踪算法：OMP.m 
function [ theta ] = CS_OMP( y,A,t )  
% 实现压缩感知OMP算法
%CS_OMP Summary of this function goes here  
%Version: 1.0 written by jbb0523 @2015-04-18  
%   Detailed explanation goes here  
%   y = Phi * x  
%   x = Psi * theta  
%   y = Phi*Psi * theta  
%   t 稀疏度
%   令 A = Phi*Psi, 则y=A*theta  
%   现在已知y和A，求theta  
%   来源：http://blog.csdn.net/jbb0523/article/details/45130793
    [y_rows,y_columns] = size(y);  
    if y_rows<y_columns  
        y = y';%y should be a column vector  
    [M,N] = size(A);%传感矩阵A为M*N矩阵  
    theta = zeros(N,1);%用来存储恢复的theta(列向量)  
    At = zeros(M,t);%用来迭代过程中存储A被选择的列  
    Pos_theta = zeros(1,t);%用来迭代过程中存储A被选择的列序号  
    r_n = y;%初始化残差(residual)为y  
    for ii=1:t%迭代t次，t为输入参数  
        product = A'*r_n;%传感矩阵A各列与残差的内积  
        [val,pos] = max(abs(product));%找到最大内积绝对值，即与残差最相关的列  
        At(:,ii) = A(:,pos);%存储这一列  
        Pos_theta(ii) = pos;%存储这一列的序号  
        A(:,pos) = zeros(M,1);%清零A的这一列，其实此行可以不要，因为它与残差正交  
        %y=At(:,1:ii)*theta，以下求theta的最小二乘解(Least Square)  
        theta_ls = (At(:,1:ii)'*At(:,1:ii))^(-1)*At(:,1:ii)'*y;%最小二乘解  
        %At(:,1:ii)*theta_ls是y在At(:,1:ii)列空间上的正交投影  
        r_n = y - At(:,1:ii)*theta_ls;%更新残差          
    theta(Pos_theta)=theta_ls;%恢复出的theta  
测试结果： 
                                    CoSaMP（Compressive Sampling Matching Pursuit）算法在MOMP算法和SP算法之间取得了平衡，它基于对原始信号的稀疏表示的假设，利用同时追踪两个系数集的方法来重构信号。具体步骤如下：首先初始化残差为原始信号，然后计算所有的内积，并选择最大的2k个系数，然后对这些系数进行逐步加和直到残差达到一定阈值或者得到了所需的系数。具体步骤如下：首先初始化残差为原始信号，然后计算所有的内积，并选择最大的k个系数，然后对这些系数进行逐步加和直到残差达到一定阈值或者得到了所需的系数。
                                    一、课题背景数据压缩技术是提高无线数据传输速度的有效措施之一。传统的数据压缩技术是基于奈奎斯特采样定律进行采样，并根据数据本身的特性降低其冗余度，从而达到压缩的目的。近年来出现的压缩感知理论（Compressed Sensing,CS）则不受制于奈奎斯特采样定律，它是采用非自适应线性投影来保持信号的原始结构，以直接采集压缩后的数据的方式，从尽量少的数据中提取尽量多的信息。本文阐述了压缩感知方法的基...
一、OMP在稀疏分解与压缩感知中的异同
1、稀疏分解要解决的问题是在冗余字典(超完备字典)A中选出k列，用这k列的线性组合近似表达待稀疏分解信号y，可以用表示为y=Aθ，求θ。
2、压缩感知重构要解决的问题是事先存在一个θ和矩阵A，然后得到y=Aθ（压缩观测），现在是在已知y和A的情况下要重构...
                                    MATLAB：matrix&laboratory两个词的组合，意为矩阵工厂(矩阵实验室)，由美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融...
                                    从上述结果可以看出，六种算法的PSNR值均在30dB以上，其中CS_IHT和CS_IRLS的PSNR值稍高，达到了33dB，而CS_GBP的PSNR值最低，仅为32.07dB。六种经典的压缩感知图像重构算法CS_CoSaMP、CS_GBP、CS_IHT、CS_IRLS、CS_OMP、CS_SP，被认为是目前最好的压缩感知图像重构算法。综上所述，通过本文的比较和分析，我们可以选择合适的压缩感知图像重构算法，以满足不同的应用需求，并提供源代码供大家学习使用。
                                    压缩感知（Compressed Sensing，CS）指出只要信号是可压缩的或在某个变换域是稀疏的,那么就可以用一个与变换基不相关的观测矩阵将变换所得高维信号投影到一个低维空间上,然后通过求解一个优化问题就可以从这些少量的投影中以高概率重构出原信号。恢复信号方式：传统采样恢复是在Nyquist采样定理的基础上，通过采样数据的sinc函数线性内插获得（在不均匀采样下则是非线性的插值恢复），而压缩感知采用的是利用信号的稀疏性，从线性观测数据中通过求解一个非线性的优化问题来恢复信号的方法。
                                    【实例简介】本工具箱包含常用的压缩感知图像重构算法，如OMP,BP,IHT,等算法，非常齐全。【实例截图】【核心代码】压缩感知图像重构算法├── BCS implemented│ ├── BCSvb.m│ ├── cameraman.pgm│ ├── demo_BCSvb_1d.m│ ├── demo_BCSvb_2d.m│ ├── demo_BCSvb_2d_small.m│...