压缩感知之常用稀疏恢复算法（贪婪算法）及matlab代码

---

前言

  压缩感知第三步是进行信号的重构，需要用到恢复重构算法。前面的文章提到过，压缩感知的恢复算法主要分为贪婪算法和凸优化算法两种，这里主要介绍贪婪算法的两种基础算法：MP算法和OMP算法及其matlab代码，并给出一些学习资料，希望可以共同进步~

---

一、匹配追踪算法(MP)

1、原理概述

  算法假定输入信号与字典库中的原子在结构上具有一定的相关性，这种相关性通过信号与原子库中原子的内积表示，即内积越大，表示信号与字典库中的这个原子的相关性越大，因此可以使用这个原子来近似表示这个信号。当然这种表示会有误差，将表示误差称为信号残差，用原信号减去这个原子，得到残差，再通过计算相关性的方式从字典库中选出一个原子表示这个残差。迭代进行上述步骤，随着迭代次数的增加，信号残差将越来越小，当满足停止条件时终止迭代，得到一组原子，及残差，将这组原子进行线性组合就能重构输入信号。

2、matlab代码

  该程序为MP算法的matlab实现，是一维信号的恢复。觉得有用就点个赞吧~

clear all
close all
%
A = [1 0.5 -1/2^0.5;
    0 (3/4)^0.5 -1/2^0.5];
y = [1,0.5]';
K = 3;
[m,n] = size(A);
 
 
% iteration
Rf(:,1) = y;
for k = 1:K
    
    for i = 1:n 
        ip(i) = abs(Rf(:,k)'*A(:,i));
    end
    j(k) = find(max(ip)==ip);
    Rf(:,k+1) = Rf(:,k) - Rf(:,k)'*A(:,j(k))*A(:,j(k));
    Rfnorm(k) = norm(Rf(:,k));
end
R = [A(:,j(1)),A(:,j(2)),A(:,j(3))];
r1 = R(:,1);
r2 = R(:,2);
r3 = R(:,3);
figure,quiver(0,0,y(1),y(2),'r');
hold,quiver(0,0,r1(1),r1(2),'b');
quiver(0,0,r2(1),r2(2),'b');
quiver(0,0,r3(1),r3(2),'b');
display(norm(Rf(:,K+1)));
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

  主要参考以下网站，并给出推荐网址：

https://www.csdn.net/tags/OtTaAg4sNjkzNTktYmxvZwO0O0OO0O0O.html
https://download.csdn.net/download/yuhongsheng1038/4888660?spm=1001.2101.3001.5697

二、正交匹配追踪算法（OMP）

1、原理概述

  正交匹配追踪算法迭代的基本思想就是每次迭代过程中从全息矩阵T中选出与测量信号s相关度（内积）最大的那一列，然后从T中去掉该列并加入到扩充矩阵Aug_t中，然后利用最小二乘法原理求出使残差r_n=s-Aug_t*aug_y最小的一个估计aug_y，然后不断的从T中减去相关列重复以上过程，直到达到迭代次数结束。参考网站如下：

https://blog.csdn.net/li200503028/article/details/14446401/
https://blog.csdn.net/scucj/article/details/7467955
https://blog.csdn.net/bingokunkun/article/details/123473709?spm=1001.2101.3001.6661.1&utm_medium=distribute.pc_relevant_t0.none-task-blog-2%7Edefault%7ECTRLIST%7ETopBlog-1.topblog&depth_1-utm_source=distribute.pc_relevant_t0.none-task-blog-2%7Edefault%7ECTRLIST%7ETopBlog-1.topblog&utm_relevant_index=1

2、matlab代码

%A-稀疏系数矩阵
%D-字典/测量矩阵（已知）
%X-测量值矩阵（已知）
%K-稀疏度
function A=OMP(D,F,X,L)
X=double(X);
[n,P]=size(X);
[n,K]=size(D);
%按列操作，分别求出每一列对应的最相关的基
for k=1:P

a=[];  %a:每一列对应的相关基的系数
x=X(:,k);   %取二维信号的每一列信号
residual=x;  %初始残差
indx=zeros(L,1);  %indx：索引集，L：测量次数
for j=1:L
	residual=double(residual);
	D=double(D);  %D转置与残差residual相乘，得到residual与每一列的内积值
	proj=D'*residual;
	pos=find(abs(proj)==max(abs(proj)));   %找内积值最大值的位置，即最相关基的position
	pos=pos(1);    %若最大值不止一个，取第一个
	indx(j)=pos;   %位置存入索引集的第j值
	a=pinv(D(:,indx(1:j)))*x;  %indx(1:j)表示第一列前j个元素;pinv：Pseudoinverse伪逆矩阵，一般用于处理长方形矩阵的求逆,得到其相关基的对  		应系数，AD=X，A=inv(D)*X;一般应该通过最小二乘来求
	residual=x-D(:,indx(1:j))*a;     %继续求残差
end

temp=zeros(K,1);   %通过上面的循环，得到第k列的相关基对应的索引位置
temp(indx)=a;
A(:,k)=temp;  %只显示非零值及其位置
end

R=A'*D';  %最终求得整个A，代入AD=X，即可求解
R=uint8(R); imshow(R);
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

  代码参考自：

https://www.cnblogs.com/AndyJee/p/3849200.html

---

总结

  本文介绍了贪婪算法中较为基础的两种算法MP和OMP算法，希望对你有所帮助~