贪婪算法中正交匹配追踪算法OMP的原理及仿真

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

前面经过几篇的基础铺垫，本篇给出正交匹配追踪(OMP)算法的MATLAB函数代码，并且给出单次测试例程代码、测量数M与重构成功概率关系曲线绘制例程代码、信号稀疏度K与重构成功概率关系曲线绘制例程代码。

0、符号说明如下：

压缩观测y=Φx，其中y为观测所得向量M×1，x为原信号N×1（M<

(1) y为观测所得向量，大小为M×1

(2)x为原信号，大小为N×1

(3)θ为K稀疏的，是信号在x在某变换域的稀疏表示

(4)Φ称为观测矩阵、测量矩阵、测量基，大小为M×N

(5)Ψ称为变换矩阵、变换基、稀疏矩阵、稀疏基、正交基字典矩阵，大小为N×N

(6)A称为测度矩阵、传感矩阵、CS信息算子，大小为M×N

上式中，一般有K<

1、OMP重构算法流程：

2、正交匹配追踪(OMP)MATLAB代码(CS_OMP.m)

[plain]view plaincopy

1.function [ theta ] = CS_OMP( y,A,t )

2.%CS_OMP Summary of this function goes here

3.%Version: 1.0 written by jbb0523 @2015-04-18

4.% Detailed explanation goes here

5.% y = Phi * x

6.% x = Psi * theta

7.% y = Phi*Psi * theta

8.% 令 A = Phi*Psi, 则y=A*theta

9.% 现在已知y和A，求theta

10. [y_rows,y_columns] = size(y);

11. if y_rows

12. y = y';%y should be a column vector

13. end

14. [M,N] = size(A);%传感矩阵A为M*N矩阵

15. theta = zeros(N,1);%用来存储恢复的theta(列向量)

16. At = zeros(M,t);%用来迭代过程中存储A被选择的列

17. Pos_theta = zeros(1,t);%用来迭代过程中存储A被选择的列序号

18. r_n = y;%初始化残差(residual)为y

19. for ii=1:t%迭代t次，t为输入参数

20. product = A'*r_n;%传感矩阵A各列与残差的内积

21. [val,pos] = max(abs(product));%找到最大内积绝对值，即与残差最相关的列

22. At(:,ii) = A(:,pos);%存储这一列

23. Pos_theta(ii) = pos;%存储这一列的序号

24. A(:,pos) = zeros(M,1);%清零A的这一列，其实此行可以不要，因为它与残差正交

25. %y=At(:,1:ii)*theta，以下求theta的最小二乘解(Least Square)

26. theta_ls = (At(:,1:ii)'*At(:,1:ii))^(-1)*At(:,1:ii)'*y;%最小二乘解

27. %At(:,1:ii)*theta_ls是y在At(:,1:ii)列空间上的正交投影

28. r_n = y - At(:,1:ii)*theta_ls;%更新残差

29. end

30. theta(Pos_theta)=theta_ls;%恢复出的theta

31.end

3、OMP单次重构测试代码(CS_Reconstuction_Test.m)

代码中，直接构造一个K稀疏的信号，所以稀疏矩阵为单位阵。

[plain]view plaincopy

1.%压缩感知重构算法测试

2.clear all;close all;clc;

3.M = 64;%观测值个数

4.N = 256;%信号x的长度

5.K = 10;%信号x的稀疏度

6.Index_K = randperm(N);

7.x = zeros(N,1);

8.x(Index_K(1:K)) = 5*randn(K,1);%x为K稀疏的，且位置是随机的

9.Psi = eye(N);%x本身是稀疏的，定义稀疏矩阵为单位阵x=Psi*theta

10.Phi = randn(M,N);%测量矩阵为高斯矩阵

11.A = Phi * Psi;%传感矩阵

12.y = Phi * x;%得到观测向量y

13.%% 恢复重构信号x

14.tic

15.theta = CS_OMP(y,A,K);

16.x_r = Psi * theta;% x=Psi * theta

17.toc

18.%% 绘图

19.figure;

20.plot(x_r,'k.-');%绘出x的恢复信号

21.hold on;

22.plot(x,'r');%绘出原信号x

23.hold off;

24.legend('Recovery','Original')

25.fprintf('\n恢复残差：');

26.norm(x_r-x)%恢复残差