贪婪算法中正交匹配追踪算法OMP的原理及仿真.docx
- 文档编号:10468387
- 上传时间:2023-02-13
- 格式:DOCX
- 页数:15
- 大小:740.11KB
贪婪算法中正交匹配追踪算法OMP的原理及仿真.docx
《贪婪算法中正交匹配追踪算法OMP的原理及仿真.docx》由会员分享,可在线阅读,更多相关《贪婪算法中正交匹配追踪算法OMP的原理及仿真.docx(15页珍藏版)》请在冰豆网上搜索。
贪婪算法中正交匹配追踪算法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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 贪婪 算法 中正 匹配 追踪 OMP 原理 仿真