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

1、压缩感知重构算法之正交匹配追踪(OMP)前面经过几篇的基础铺垫，本篇给出正交匹配追踪(OMP)算法的MATLAB函数代码，并且给出单次测试例程代码、测量数M与重构成功概率关系曲线绘制例程代码、信号稀疏度K与重构成功概率关系曲线绘制例程代码。0、符号说明如下： 压缩观测y=x，其中y为观测所得向量M1，x为原信号N1（MN）。x一般不是稀疏的，但在某个变换域是稀疏的，即x=，其中为K稀疏的，即只有K个非零项。此时y=，令A=，则y=A。 (1)y为观测所得向量，大小为M1 (2)x为原信号，大小为N1 (3)为K稀疏的，是信号在x在某变换域的稀疏表示 (4)称为观测矩阵、测量矩阵、测量基，大小为

2、MN (5)称为变换矩阵、变换基、稀疏矩阵、稀疏基、正交基字典矩阵，大小为NN (6)A称为测度矩阵、传感矩阵、CS信息算子，大小为MN上式中，一般有KMN，后面三个矩阵各个文献的叫法不一，以后我将称为测量矩阵、将称为稀疏矩阵、将A称为传感矩阵。1、OMP重构算法流程：2、正交匹配追踪(OMP)MATLAB代码(CS_OMP.m)plainview plaincopy1. functiontheta=CS_OMP(y,A,t)2. %CS_OMPSummaryofthisfunctiongoeshere3. %Version:1.0writtenbyjbb05232015-04-184. %D

3、etailedexplanationgoeshere5. %y=Phi*x6. %x=Psi*theta7. %y=Phi*Psi*theta8. %令A=Phi*Psi,则y=A*theta9. %现在已知y和A，求theta10. y_rows,y_columns=size(y);11. ify_rowsy_columns12. y=y;%yshouldbeacolumnvector13. end14. M,N=size(A);%传感矩阵A为M*N矩阵15. theta=zeros(N,1);%用来存储恢复的theta(列向量)16. At=zeros(M,t);%用来迭代过程中存储A被选

4、择的列17. Pos_theta=zeros(1,t);%用来迭代过程中存储A被选择的列序号18. r_n=y;%初始化残差(residual)为y19. forii=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(:

5、,1:ii)*theta，以下求theta的最小二乘解(LeastSquare)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. end30. theta(Pos_theta)=theta_ls;%恢复出的theta31. end3、OMP单次重构测试代码(CS_Reconstuction_Test.m) 代码中，直接构造一个K稀疏的信号，所以稀疏矩阵为单位阵

6、。plainview plaincopy1. %压缩感知重构算法测试2. clearall;closeall;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*theta10. Phi=randn(M,N);%测量矩阵为高斯矩阵11. A=Phi*Psi;%传感矩阵12. y=Phi*x;%得到观测向

7、量y13. %恢复重构信号x14. tic15. theta=CS_OMP(y,A,K);16. x_r=Psi*theta;%x=Psi*theta17. toc18. %绘图19. figure;20. plot(x_r,k.-);%绘出x的恢复信号21. holdon;22. plot(x,r);%绘出原信号x23. holdoff;24. legend(Recovery,Original)25. fprintf(n恢复残差：);26. norm(x_r-x)%恢复残差运行结果如下：（信号为随机生成，所以每次结果均不一样）1）图：2）Command WindowsElapsed time

8、 is 0.849710 seconds.恢复残差：ans = 5.5020e-0154、测量数M与重构成功概率关系曲线绘制例程代码plainview plaincopy1. %压缩感知重构算法测试CS_Reconstuction_MtoPercentage.m2. %绘制参考文献中的Fig.13. %参考文献：JoelA.TroppandAnnaC.Gilbert4. %SignalRecoveryFromRandomMeasurementsViaOrthogonalMatching5. %Pursuit，IEEETRANSACTIONSONINFORMATIONTHEORY,VOL.53,

9、NO.12,6. %DECEMBER2007.7. %Elapsedtimeis1171.606254seconds.(20150418night)8. clearall;closeall;clc;9. %参数配置初始化10. CNT=1000;%对于每组(K,M,N)，重复迭代次数11. N=256;%信号x的长度12. Psi=eye(N);%x本身是稀疏的，定义稀疏矩阵为单位阵x=Psi*theta13. K_set=4,12,20,28,36;%信号x的稀疏度集合14. Percentage=zeros(length(K_set),N);%存储恢复成功概率15. %主循环，遍历每组(K

10、,M,N)16. tic17. forkk=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. formm=1:length(M_set)22. M=M_set(mm);%本次观测值个数23. P=0;24. forcnt=1:CNT%每个观测值个数均运行CNT次25. Index_K=randperm(N);26. x=zeros(N,1);27. x(Index_K(1:K)

11、=5*randn(K,1);%x为K稀疏的，且位置是随机的28. Phi=randn(M,N);%测量矩阵为高斯矩阵29. A=Phi*Psi;%传感矩阵30. y=Phi*x;%得到观测向量y31. theta=CS_OMP(y,A,K);%恢复重构信号theta32. x_r=Psi*theta;%x=Psi*theta33. ifnorm(x_r-x)1e-6%如果残差小于1e-6则认为恢复成功34. P=P+1;35. end36. end37. PercentageK(mm)=P/CNT*100;%计算恢复概率38. end39. Percentage(kk,1:length(M_s

12、et)=PercentageK;40. end41. toc42. saveMtoPercentage1000%运行一次不容易，把变量全部存储下来43. %绘图44. S=-ks;-ko;-kd;-kv;-k*;45. figure;46. forkk=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. holdon;52. end53. holdoff;54. xlim(0256);

13、55. legend(K=4,K=12,K=20,K=28,K=36);56. xlabel(Numberofmeasurements(M);57. ylabel(Percentagerecovered);58. title(Percentageofinputsignalsrecoveredcorrectly(N=256)(Gaussian); 本程序在联想ThinkPadE430C笔记本（4GB DDR3内存，i5-3210）上运行共耗时1171.606254秒，程序中将所有数据均通过“save MtoPercentage1000”存储了下来，以后可以再对数据进行分析，只需“load Mto

14、Percentage1000”即可。 程序运行结果比文献1的Fig.1要好，原因不详。本程序运行结果：文献1中的Fig.1：5、信号稀疏度K与重构成功概率关系曲线绘制例程代码plainview plaincopy1. %压缩感知重构算法测试CS_Reconstuction_KtoPercentage.m2. %绘制参考文献中的Fig.23. %参考文献：JoelA.TroppandAnnaC.Gilbert4. %SignalRecoveryFromRandomMeasurementsViaOrthogonalMatching5. %Pursuit，IEEETRANSACTIONSONINFO

15、RMATIONTHEORY,VOL.53,NO.12,6. %DECEMBER2007.7. %Elapsedtimeis1448.966882seconds.(20150418night)8. clearall;closeall;clc;9. %参数配置初始化10. CNT=1000;%对于每组(K,M,N)，重复迭代次数11. N=256;%信号x的长度12. Psi=eye(N);%x本身是稀疏的，定义稀疏矩阵为单位阵x=Psi*theta13. M_set=52,100,148,196,244;%测量值集合14. Percentage=zeros(length(M_set),N);%存

16、储恢复成功概率15. %主循环，遍历每组(K,M,N)16. tic17. formm=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. forkk=1:length(K_set)22. K=K_set(kk);%本次信号x的稀疏度K23. P=0;24. forcnt=1:CNT%每个观测值个数均运行CNT次25. Index_K=randp

17、erm(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;%得到观测向量y31. theta=CS_OMP(y,A,K);%恢复重构信号theta32. x_r=Psi*theta;%x=Psi*theta33. ifnorm(x_r-x)1e-6%如果残差小于1e-6则认为恢复成功34. P=P+1;35. end36. end37. PercentageM(kk)=P/CNT*100;%计

18、算恢复概率38. end39. Percentage(mm,1:length(K_set)=PercentageM;40. end41. toc42. saveKtoPercentage1000test%运行一次不容易，把变量全部存储下来43. %绘图44. S=-ks;-ko;-kd;-kv;-k*;45. figure;46. formm=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. holdon;52. end53. holdoff;54. xlim(0125);55. legend(M=52,M=100,M=148,M=196,M=244);56. xlabel(Sparsitylevel(K);57. ylabel(Percentagerecovered);58. title(Percentageofinputsignalsrecoveredcorrectly(N=256)(Gaussian); 本程序在联想ThinkPadE430C笔记本（4GB DDR3内存，i5-3210）上运行共耗时1448.966882秒，程序中将所有数