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

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

《贪婪算法中正交匹配追踪算法OMP的原理及仿真.docx》由会员分享，可在线阅读，更多相关《贪婪算法中正交匹配追踪算法OMP的原理及仿真.docx（13页珍藏版）》请在冰豆网上搜索。

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

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

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

0、符号说明如下：

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

x一般不是稀疏的，但在某个变换域Ψ是稀疏的，即x=Ψθ，其中θ为K稀疏的，即θ只有K个非零项。

此时y=ΦΨθ，令A=ΦΨ，则y=Aθ。

（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] viewplaincopy

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] viewplaincopy

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）%恢复残差  

运行结果如下：

（信号为随机生成，所以每次结果均不一样）

1）图：

2）CommandWindows

Elapsedtimeis0.849710seconds.

恢复残差：

ans=

 5.5020e-015

4、测量数M与重构成功概率关系曲线绘制例程代码

[plain] viewplaincopy

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

2.%   绘制参考文献中的Fig.1  

3.%   参考文献：

Joel A. Tropp and Anna C. Gilbert   

4.%   Signal Recovery From Random Measurements Via Orthogonal Matching  

5.%   Pursuit，IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 53, NO. 12,  

6.%   DECEMBER 2007.  

7.%   Elapsed time is 1171.606254 seconds.（@20150418night）  

8.clear all;close all;clc;  

9.%% 参数配置初始化  

10.CNT = 1000;%对于每组（K,M,N），重复迭代次数  

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

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

13.K_set = [4,12,20,28,36];%信号x的稀疏度集合  

14.Percentage = zeros（length（K_set）,N）;%存储恢复成功概率  

15.%% 主循环，遍历每组（K,M,N）  

16.tic  

17.for kk = 1:

length（K_set）  

18.    K = K_set（kk）;%本次稀疏度  

19.    M_set = K:

5:

N;%M没必要全部遍历，每隔5测试一个就可以了  

20.    PercentageK = zeros（1,length（M_set））;%存储此稀疏度K下不同M的恢复成功概率  

21.    for mm = 1:

length（M_set）  

22.       M = M_set（mm）;%本次观测值个数  

23.       P = 0;  

24.       for cnt = 1:

CNT %每个观测值个数均运行CNT次  

25.            Index_K = randperm（N）;  

26.            x = zeros（N,1）;  

27.            x（Index_K（1:

K）） = 5*randn（K,1）;%x为K稀疏的，且位置是随机的                  

28.            Phi = randn（M,N）;%测量矩阵为高斯矩阵  

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

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

31.            theta = CS_OMP（y,A,K）;%恢复重构信号theta  

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

33.            if norm（x_r-x）<1e-6%如果残差小于1e-6则认为恢复成功  

34.                P = P + 1;  

35.            end  

36.       end  

37.       PercentageK（mm） = P/CNT*100;%计算恢复概率  

38.    end  

39.    Percentage（kk,1:

length（M_set）） = PercentageK;  

40.end  

41.toc  

42.save MtoPercentage1000 %运行一次不容易，把变量全部存储下来  

43.%% 绘图  

44.S = ['-ks';'-ko';'-kd';'-kv';'-k*'];  

45.figure;  

46.for kk = 1:

length（K_set）  

47.    K = K_set（kk）;  

48.    M_set = K:

5:

N;  

49.    L_Mset = length（M_set）;  

50.    plot（M_set,Percentage（kk,1:

L_Mset）,S（kk,:

））;%绘出x的恢复信号  

51.    hold on;  

52.end  

53.hold off;  

54.xlim（[0 256]）;  

55.legend（'K=4','K=12','K=20','K=28','K=36'）;  

56.xlabel（'Number of measurements（M）'）;  

57.ylabel（'Percentage recovered'）;  

58.title（'Percentage of input signals recovered correctly（N=256）（Gaussian）'）;  

    本程序在联想ThinkPadE430C笔记本（4GBDDR3内存，i5-3210）上运行共耗时1171.606254秒，程序中将所有数据均通过“saveMtoPercentage1000”存储了下来，以后可以再对数据进行分析，只需“loadMtoPercentage1000”即可。

    程序运行结果比文献[1]的Fig.1要好，原因不详。

本程序运行结果：

文献[1]中的Fig.1：

5、信号稀疏度K与重构成功概率关系曲线绘制例程代码

[plain] viewplaincopy

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

2.%   绘制参考文献中的Fig.2  

3.%   参考文献：

Joel A. Tropp and Anna C. Gilbert   

4.%   Signal Recovery From Random Measurements Via Orthogonal Matching  

5.%   Pursuit，IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 53, NO. 12,  

6.%   DECEMBER 2007.  

7.%   Elapsed time is 1448.966882 seconds.（@20150418night）  

8.clear all;close all;clc;  

9.%% 参数配置初始化  

10.CNT = 1000;%对于每组（K,M,N），重复迭代次数  

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

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

13.M_set = [52,100,148,196,244];%测量值集合  

14.Percentage = zeros（length（M_set）,N）;%存储恢复成功概率  

15.%% 主循环，遍历每组（K,M,N）  

16.tic  

17.for mm = 1:

length（M_set）  

18.    M = M_set（mm）;%本次测量值个数  

19.    K_set = 1:

5:

ceil（M/2）;%信号x的稀疏度K没必要全部遍历，每隔5测试一个就可以了  

20.    PercentageM = zeros（1,length（K_set））;%存储此测量值M下不同K的恢复成功概率  

21.    for kk = 1:

length（K_set）  

22.       K = K_set（kk）;%本次信号x的稀疏度K  

23.       P = 0;  

24.       for cnt = 1:

CNT %每个观测值个数均运行CNT次  

25.            Index_K = randperm（N）;  

26.            x = zeros（N,1）;  

27.            x（Index_K（1:

K）） = 5*randn（K,1）;%x为K稀疏的，且位置是随机的                  

28.            Phi = randn（M,N）;%测量矩阵为高斯矩阵  

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

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

31.            theta = CS_OMP（y,A,K）;%恢复重构信号theta  

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

33.            if norm（x_r-x）<1e-6%如果残差小于1e-6则认为恢复成功  

34.                P = P + 1;  

35.            end  

36.       end  

37.       PercentageM（kk） = P/CNT*100;%计算恢复概率  

38.    end  

39.    Percentage（mm,1:

length（K_set）） = PercentageM;  

40.end  

41.toc  

42.save KtoPercentage1000test %运行一次不容易，把变量全部存储下来  

43.%% 绘图  

44.S = ['-ks';'-ko';'-kd';'-kv';'-k*'];  

45.figure;  

46.for mm = 1:

length（M_set）  

47.    M = M_set（mm）;  

48.    K_set = 1:

5:

ceil（M/2）;  

49.    L_Kset = length（K_set）;  

50.    plot（K_set,Percentage（mm,1:

L_Kset）,S（mm,:

））;%绘出x的恢复信号  

51.    hold on;  

52.end  

53.hold off;  

54.xlim（[0 125]）;  

55.legend（'M=52','M=100','M=148','M=196','M=244'）;  

56.xlabel（'Sparsity level（K）'）;  

57.ylabel（'Percentage recovered'）;  

58.title（'Percentage of input signals recovered correctly（N=256）（Gaussian）'）;  

    本程序在联想ThinkPadE430C笔记本（4GBDDR3内存，i5-3210）上运行共耗时1448.966882秒，程序中将所有数据均通过“saveKtoPercentage1000”存储了下来，以后可以再对数据进行分析，只需“loadKtoPercentage1000”即可。

    程序运行结果比文献[1]的Fig.2要好，原因不详。

本程序运行结果：

文献[1]中的Fig.2：

6、参考文献

【1】JoelA.TroppandAnnaC.Gilbert.SignalRecoveryFromRandomMeasurementsViaOrthogonalMatchingPursuit[J].IEEETransactionsonInformationTheory,VOL.53,NO.12,DECEMBER2007.

【2】Y.C.Pati,R.Rezaiifar,andP.S.Krishnaprasad.OrthogonalMatchingPursuit-RecursiveFunctionApproximationwithApplicationstowaveletdecomposition,Proc.27thAnnu.AsilomarConf.Signals,Systems,andComputers,PacificGrove,CA,Nov.1993,vol.1,pp40-44.

【3】沙威.CS_OMP，http:

//www.eee.hku.hk/~wsha/Freecode/Files/CS_OMP.zip