OMP算法代码学习

article/2025/4/30 15:36:48
      版权声明:本文为博主原创文章,转载请注明出处,谢谢!          https://blog.csdn.net/jbb0523/article/details/45130793        </div><link rel="stylesheet" href="https://csdnimg.cn/release/phoenix/template/css/ck_htmledit_views-f57960eb32.css"><link rel="stylesheet" href="https://csdnimg.cn/release/phoenix/template/css/ck_htmledit_views-f57960eb32.css"><div class="htmledit_views" id="content_views">

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

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

0、符号说明如下:

        压缩观测y=Φx,其中y为观测所得向量M×1,x为原信号N×1(M<<N)。x一般不是稀疏的,但在某个变换域Ψ是稀疏的,即x=Ψθ,其中θ为K稀疏的,即θ只有K个非零项。此时y=ΦΨθ,令A=ΦΨ,则y=

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

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

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

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

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

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

上式中,一般有K<<M<<N,后面三个矩阵各个文献的叫法不一,以后我将Φ称为测量矩阵、将Ψ称为稀疏矩阵、将A称为传感矩阵

1、OMP重构算法流程:



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

function [ theta ] = CS_OMP( y,A,t )
%CS_OMP Summary of this function goes here
%Version: 1.0 written by jbb0523 @2015-04-18
%   Detailed explanation goes here
%   y = Phi * x
%   x = Psi * theta
%	y = Phi*Psi * theta
%   令 A = Phi*Psi, 则y=A*theta
%   现在已知y和A,求theta[y_rows,y_columns] = size(y);if y_rows<y_columnsy = y';%y should be a column vectorend[M,N] = size(A);%传感矩阵A为M*N矩阵theta = zeros(N,1);%用来存储恢复的theta(列向量)At = zeros(M,t);%用来迭代过程中存储A被选择的列Pos_theta = zeros(1,t);%用来迭代过程中存储A被选择的列序号r_n = y;%初始化残差(residual)为yfor ii=1:t%迭代t次,t为输入参数product = A'*r_n;%传感矩阵A各列与残差的内积[val,pos] = max(abs(product));%找到最大内积绝对值,即与残差最相关的列At(:,ii) = A(:,pos);%存储这一列Pos_theta(ii) = pos;%存储这一列的序号A(:,pos) = zeros(M,1);%清零A的这一列,其实此行可以不要,因为它与残差正交%y=At(:,1:ii)*theta,以下求theta的最小二乘解(Least Square)theta_ls = (At(:,1:ii)'*At(:,1:ii))^(-1)*At(:,1:ii)'*y;%最小二乘解%At(:,1:ii)*theta_lsyAt(:,1:ii)列空间上的正交投影r_n = y - At(:,1:ii)*theta_ls;%更新残差        endtheta(Pos_theta)=theta_ls;%恢复出的theta
end


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

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

%压缩感知重构算法测试
clear all;close all;clc;
M = 64;%观测值个数
N = 256;%信号x的长度
K = 10;%信号x的稀疏度
Index_K = randperm(N);
x = zeros(N,1);
x(Index_K(1:K)) = 5*randn(K,1);%x为K稀疏的,且位置是随机的
Psi = eye(N);%x本身是稀疏的,定义稀疏矩阵为单位阵x=Psi*theta
Phi = randn(M,N);%测量矩阵为高斯矩阵
A = Phi * Psi;%传感矩阵
y = Phi * x;%得到观测向量y
%% 恢复重构信号x
tic
theta = CS_OMP(y,A,K);
x_r = Psi * theta;% x=Psi * theta
toc
%% 绘图
figure;
plot(x_r,'k.-');%绘出x的恢复信号
hold on;
plot(x,'r');%绘出原信号x
hold off;
legend('Recovery','Original')
fprintf('\n恢复残差:');
norm(x_r-x)%恢复残差

运行结果如下:(信号为随机生成,所以每次结果均不一样)

1)图:


2)Command Windows

Elapsed time is 0.849710 seconds.
恢复残差:
ans =
  5.5020e-015

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

%压缩感知重构算法测试CS_Reconstuction_MtoPercentage.m
%   绘制参考文献中的Fig.1
%   参考文献:Joel A. Tropp and Anna C. Gilbert 
%   Signal Recovery From Random Measurements Via Orthogonal Matching
%   Pursuit,IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 53, NO. 12,
%   DECEMBER 2007.
%   Elapsed time is 1171.606254 seconds.(@20150418night)
clear all;close all;clc;
%% 参数配置初始化
CNT = 1000;%对于每组(K,M,N),重复迭代次数
N = 256;%信号x的长度
Psi = eye(N);%x本身是稀疏的,定义稀疏矩阵为单位阵x=Psi*theta
K_set = [4,12,20,28,36];%信号x的稀疏度集合
Percentage = zeros(length(K_set),N);%存储恢复成功概率
%% 主循环,遍历每组(K,M,N)
tic
for kk = 1:length(K_set)K = K_set(kk);%本次稀疏度M_set = K:5:N;%M没必要全部遍历,每隔5测试一个就可以了PercentageK = zeros(1,length(M_set));%存储此稀疏度K下不同M的恢复成功概率for mm = 1:length(M_set)M = M_set(mm);%本次观测值个数P = 0;for cnt = 1:CNT %每个观测值个数均运行CNT次Index_K = randperm(N);x = zeros(N,1);x(Index_K(1:K)) = 5*randn(K,1);%x为K稀疏的,且位置是随机的                Phi = randn(M,N);%测量矩阵为高斯矩阵A = Phi * Psi;%传感矩阵y = Phi * x;%得到观测向量ytheta = CS_OMP(y,A,K);%恢复重构信号thetax_r = Psi * theta;% x=Psi * thetaif norm(x_r-x)<1e-6%如果残差小于1e-6则认为恢复成功P = P + 1;endendPercentageK(mm) = P/CNT*100;%计算恢复概率endPercentage(kk,1:length(M_set)) = PercentageK;
end
toc
save MtoPercentage1000 %运行一次不容易,把变量全部存储下来
%% 绘图
S = ['-ks';'-ko';'-kd';'-kv';'-k*'];
figure;
for kk = 1:length(K_set)K = K_set(kk);M_set = K:5:N;L_Mset = length(M_set);plot(M_set,Percentage(kk,1:L_Mset),S(kk,:));%绘出x的恢复信号hold on;
end
hold off;
xlim([0 256]);
legend('K=4','K=12','K=20','K=28','K=36');
xlabel('Number of measurements(M)');
ylabel('Percentage recovered');
title('Percentage of input signals recovered correctly(N=256)(Gaussian)');

        本程序在联想ThinkPadE430C笔记本(4GB DDR3内存,i5-3210)上运行共耗时1171.606254秒,程序中将所有数据均通过“save MtoPercentage1000”存储了下来,以后可以再对数据进行分析,只需“load MtoPercentage1000”即可。

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

本程序运行结果:


文献[1]中的Fig.1:


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

%压缩感知重构算法测试CS_Reconstuction_KtoPercentage.m
%   绘制参考文献中的Fig.2
%   参考文献:Joel A. Tropp and Anna C. Gilbert 
%   Signal Recovery From Random Measurements Via Orthogonal Matching
%   Pursuit,IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 53, NO. 12,
%   DECEMBER 2007.
%   Elapsed time is 1448.966882 seconds.(@20150418night)
clear all;close all;clc;
%% 参数配置初始化
CNT = 1000;%对于每组(K,M,N),重复迭代次数
N = 256;%信号x的长度
Psi = eye(N);%x本身是稀疏的,定义稀疏矩阵为单位阵x=Psi*theta
M_set = [52,100,148,196,244];%测量值集合
Percentage = zeros(length(M_set),N);%存储恢复成功概率
%% 主循环,遍历每组(K,M,N)
tic
for mm = 1:length(M_set)M = M_set(mm);%本次测量值个数K_set = 1:5:ceil(M/2);%信号x的稀疏度K没必要全部遍历,每隔5测试一个就可以了PercentageM = zeros(1,length(K_set));%存储此测量值M下不同K的恢复成功概率for kk = 1:length(K_set)K = K_set(kk);%本次信号x的稀疏度KP = 0;for cnt = 1:CNT %每个观测值个数均运行CNT次Index_K = randperm(N);x = zeros(N,1);x(Index_K(1:K)) = 5*randn(K,1);%x为K稀疏的,且位置是随机的                Phi = randn(M,N);%测量矩阵为高斯矩阵A = Phi * Psi;%传感矩阵y = Phi * x;%得到观测向量ytheta = CS_OMP(y,A,K);%恢复重构信号thetax_r = Psi * theta;% x=Psi * thetaif norm(x_r-x)<1e-6%如果残差小于1e-6则认为恢复成功P = P + 1;endendPercentageM(kk) = P/CNT*100;%计算恢复概率endPercentage(mm,1:length(K_set)) = PercentageM;
end
toc
save KtoPercentage1000test %运行一次不容易,把变量全部存储下来
%% 绘图
S = ['-ks';'-ko';'-kd';'-kv';'-k*'];
figure;
for mm = 1:length(M_set)M = M_set(mm);K_set = 1:5:ceil(M/2);L_Kset = length(K_set);plot(K_set,Percentage(mm,1:L_Kset),S(mm,:));%绘出x的恢复信号hold on;
end
hold off;
xlim([0 125]);
legend('M=52','M=100','M=148','M=196','M=244');
xlabel('Sparsity level(K)');
ylabel('Percentage recovered');
title('Percentage of input signals recovered correctly(N=256)(Gaussian)');

        本程序在联想ThinkPadE430C笔记本(4GB DDR3内存,i5-3210)上运行共耗时1448.966882秒,程序中将所有数据均通过“save KtoPercentage1000”存储了下来,以后可以再对数据进行分析,只需“load KtoPercentage1000”即可。

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

本程序运行结果:


文献[1]中的Fig.2:


6、参考文献

【1】Joel A. Tropp and Anna C. Gilbert. Signal Recovery From Random Measurements Via Orthogonal Matching Pursuit[J]. IEEETransactions on Information Theory, VOL. 53, NO. 12, DECEMBER 2007.

【2】Y.C.Pati, R.Rezaiifar,and P.S.Krishnaprasad. Orthogonal Matching Pursuit-Recursive FunctionApproximation with Applications to wavelet decomposition, Proc. 27thAnnu. Asilomar Conf. Signals, Systems, and Computers, Pacific Grove, CA, Nov.1993,vol.1,pp40-44.

【3】沙威.CS_OMP,http://www.eee.hku.hk/~wsha/Freecode/Files/CS_OMP.zip


http://chatgpt.dhexx.cn/article/sAKXCZFY.shtml

相关文章

美团饿了么返利公众号小程序搭建(付源码)

外卖返利小程序源码; 轻松部署搭建&#xff0c;小程序服务号数据互通&#xff1b; 对接美团官方; 佣金比例自定义分配; 三级分佣&#xff0c;所有资金数据一目了然&#xff1b; 拉新立减最低4.9元购月卡&#xff1b; 签到20天免费领取会员卡&#xff1b; 提现秒到账&#xff01…

OMP算法笔记

OMP算法笔记 OMP算法整理&#xff08;以备自己后期查阅&#xff0c;集合了几篇博主的文章&#xff09; &#xff08;1&#xff09;数理知识基础–投影矩阵 详见&#xff1a; 作者&#xff1a;nineheaded_bird 来源&#xff1a;CSDN 原文&#xff1a;https://blog.csdn.net/teng…

使用python实现微信小程序自动签到2.0

微信小程序自动签到 功能描述目标输出包管理 程序的结构设计步骤1步骤2步骤3步骤4 代码实现使用findler抓包工具查看请求类型再次使用findler抓包&#xff0c;查看请求内容使用多线程完成多用户提交的功能使用itchat第三方库实现微信自动回复 将程序部署到服务器中使用scp命令将…

并行程序设计——OMP编程

并行程序设计——OMP编程 实验一 实验内容 分别实现课件中的梯形积分法的Pthread、OpenMP版本&#xff0c;熟悉并掌握OpenMP编程方法&#xff0c;探讨两种编程方式的异同。 实验代码 OpenMP编程 #include <stdio.h> #include <stdlib.h> #include <omp.h&g…

利用中国知网快速自动生成参考文献

1.打开知网 2.输入引用的文献名称 3.点击文章 4.点击“导出/参考文献” 5.自动生成参考文献

知网导出外文参考文献格式和下载文章(2019.5)

如果不弹出页面&#xff0c;是网络的原因&#xff0c;等一等

谷歌学术、中国知网生成参考文献

谷歌学术 step1.登陆谷歌学术 step2.查找需要的文献 step3.点击引用标志 step4.生成相关引用 step5.选择不同的标准复制粘贴 中国知网 step1.在知网搜索需要的论文 step2.点击导出参考文献 step3.生成参考文献引用

中国知网文献引用导入EndNote9.X,Web of science导入endnote以及谷歌学术导入endnote图文详解,全网最细版本适用EndNote9.x,Endnote20版本

文章目录 一、EndNote导入文献的以下几种格式1.1 中国知网1.2web of science1.3 谷歌学术 一、EndNote导入文献的以下几种格式 all as we konow&#xff0c;引用参考文献也就是如下三个&#xff0c;那么分别导入endnote改怎么使用呢&#xff1f; 1.1 中国知网 在你想要的文献里…

在CNKI上导出TXT文件

在使用CiteSpace之前要先下载数据源&#xff0c;今天就来讲一讲从CNKI上导出txt文件。 1、从学校官网进入中国知网CNKI&#xff0c;单击高级检索 2、输入关键字&#xff0c;可以选择组合输入&#xff0c;单击搜索 3、在每页显示处选择50 4、勾选所有所有记录&#xff08;每次导…

EndNote20导入知网文献和导出BibTeX于TexStudio

第一步&#xff1a;在知网官网中找到一篇论文并点击引号键 第二步&#xff1a;点击EndNote下载该txt文件 第三步&#xff1a;在EndNote20中点击File->Import->File&#xff0c;引入刚刚下载的txt文件 注意Import Option选择的是EndNote Import&#xff0c;点击import导入…

知网导出之Excel

背景&#xff1a;需要整理知网XX方向的论文&#xff0c;包括题目&#xff0c;内容&#xff0c;发表时间等等信息。但是一个个写也太麻烦了&#xff0c;所以&#xff0c;一些不需要人总结的东西&#xff0c;就直接导出好了。 流程&#xff1a; 选择好文献之后&#xff0c;选择自…

知网如何快速引用参考文献

1.在知网搜索界面&#xff0c;搜索自己想搜的内容&#xff0c;然后点击下图引号位置 2.在弹框中选择你要引用的格式 复制粘贴到自己论文中即可

知网 BibTeX自动生成(使用BibTeX引用中文参考文献)

前言 谷歌学术具备生成英文文献的bibtex文献引用代码的功能&#xff0c;而知网里不具备生成中文文献的bibtex引用代码的功能。因此&#xff0c;本文将生成中文文献bibtex引用代码的操作过程简单记录便于自己再次翻阅&#xff0c;操作方法源自知乎作者 上官无忌 &#xff0c;具…

知网导出引用文件,插入到Endnote管理文献

直接选endnote是txt格式 改格式也没用&#xff0c;endnote改也不行 这里选择refworks 导出txt之后 选择refworks import即可

Web of Science如何导出参考文献

Web of Science&#xff08;WOS&#xff09;是大型综合性、多学科、核心期刊引文索引数据库。 Endnote是一款可以有效管理参考文献的软件&#xff0c;并且在论文写作时&#xff0c;利用endnote可以很方便的修改引用格式。 1.进入论文界面 2.点击查看PDF,然后点击红圈 3. 可以…

Endnote导出目标期刊的参考文献的格式

Endnote导出目标期刊的参考文献的格式 1、知网文献与endnote的操作流程2、其他的参考文件来源&#xff08;如bib、ris等&#xff09;及导出格式选择 1、知网文献与endnote的操作流程 导出的格式为txt格式&#xff0c;这时候打开endnote&#xff0c;导入choose你的txt文件就可以…

知网论文参考文献导入到Endnote方法

第一步下载参考文献文件第二步&#xff1a;右键->打开方式->选择endnote 3. 导入结果

知网获取论文参考文献

知网获取论文参考文献 进入知网搜索相应材料普通检索高级检索 选择相应的文献点击右上角左边双引号“凑”参考文献 进入知网 中国知网官方网址&#xff1a;https://www.cnki.net/ 搜索相应材料 搜素一般可分为普通检索和高级检索。 一般而言&#xff0c;普通检索即可完成我…

Endnote 导出英文、中文(知网)参考文献进入Word

1、英文文献 从Google Scholar 搜索需要的参考文献&#xff0c;然后点击“引用”按钮&#xff0c;导出Endnote的格式&#xff0c;例如scholar.enw。 在Endnote中File-->Import-->File...-->Import File-->Import 参考文献导入完毕 进行参考文献在word中的导出…

将CNKI知网文献条目导出,并导入到Endnote内

将CNKI知网文献条目导出&#xff0c;并导入到Endnote内 目录 将CNKI知网文献条目导出&#xff0c;并导入到Endnote内一、从知网上导出参考文献二、将知网导出的参考文献导入到Endnote 一、从知网上导出参考文献 从知网上导出参考文献过程和步骤如图1所示。 图1 导出的参考文献…