压缩感知重构算法之正交匹配追踪(omp)及其matlab实现

2023-10-07 18:10

本文主要是介绍压缩感知重构算法之正交匹配追踪(omp)及其matlab实现,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

压缩感知之OMP恢复算法

1、基本思想

  y=Φx   x=Ψθ
  正交匹配追踪算法的本质思想是,以贪婪迭代的方式选择测量矩阵Φ的列,使得在每次迭代中所选择的列与当前的冗余向量最大程度地相关,从测量向量中减去相关部分并反复迭代,直到迭代次数达到稀疏度K,强制迭代停止。

2、算法步骤

  输入:(1)M*N的感知矩阵A,其中M远远小于N,A=Φ*Ψ。
     (2)长度为M的数据向量b,即测量值y。
  输出:长度为N的重建向量 xˆ ,满足y=Ax。
  初始化:残差r0=y,重建信号x0=0,索引集Λ0=Φ,迭代次数n=2*K,计数器k=0。
  步骤1:计算残差和感知矩阵A的每一列的投影系数(内积值) ck=ATrk1
  步骤2:找出ck中元素最大的元素 ck=max{ck} 以及对应的位置pos;
  步骤3:更新索引集 Λk=Λk1{pos}, 以及原子集合 AΛK=AΛk1{A(:pos)}
  步骤4:利用最小二乘求得近似解 xk=(ATΛkAΛk1)1ATΛky
  步骤5:更新余量 rk=yAxk
  步骤6:判断迭代是否满足停止条件,满足则停止 xˆ=xk,r=rk , 输出 xˆ,r ,否则转步骤1。

3、仿真验证

  3.1 一维时间稀疏信号
  首先进行一维时间稀疏信号的恢复,信号长度为512,稀疏度选取10、20、30、40、50,matlab代码如下:

clc;clear;close all
%%  1. 时域测试信号生成
CNT = 100;                                       %对于每组(K,M,N),重复迭代次数  
N=512;                                           %信号长度
K_set= [10,20,30,40,50];                             %信号x的稀疏度集合 
Percentage = zeros(length(K_set),N);             %存储恢复成功概率
for kk=1:length(K_set)K=K_set(kk);                                      %本次稀疏度M_set=1:5:N;                                          %测量数每隔五个取一次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);                              %将1-N随机打乱 行向量x=zeros(N,1);x(Index_K(1:K))=5*randn(K,1);                     %x为K稀疏的,且位置是随机的%%  2.  时域信号压缩传感Phi=randn(M,N);                                   %  测量矩阵(高斯分布白噪声)Phi = Phi./repmat(sqrt(sum(Phi.^2,1)),[floor(M),1]); %正则化y=Phi*x;                                        %  获得线性测量 %%  3.  正交匹配追踪法重构信号(本质上是L_1范数最优化问题)
n=2*K;                                            %  算法迭代次数(m>=K)
Psi=eye(N);                                       %  单位矩阵为正变换矩阵
A=Phi*Psi;                                        %  恢复矩阵(测量矩阵*正交反变换矩阵)hat_y=zeros(1,N);                                 %  待重构的变换域里的向量                     
Aug_t=[];                                         %  增量矩阵(初始值为空矩阵)
r0=y;                                             %  残差值for times=1:n;                                    %  迭代次数(有噪声的情况下,该迭代次数为K)for col=1:N;                                  %  恢复矩阵的所有列向量    步骤1product(col)=abs(A(:,col)'*r0);           %  恢复矩阵的列向量和残差的投影系数(内积值) end[val,pos]=max(product);                       %  最大投影系数对应的位置  步骤2Aug_t=[Aug_t,A(:,pos)];                       %  矩阵扩充               步骤3 更新原子集合A(:,pos)=zeros(M,1);                          %  选中的列置零(实质上应该去掉,为了简单我把它置零)aug_y=(Aug_t'*Aug_t)^(-1)*Aug_t'*y;           %  最小二乘,使残差最小     步骤4 求近似解r0=y-Aug_t*aug_y;                             %  残差                   步骤5 更新余量pos_array(times)=pos;                         %  纪录最大投影系数的位置  步骤3 更新索引集
end
hat_y(pos_array)=aug_y;                           %  重构的变换域里的向量
hat_x=real(Psi'*hat_y.');                         %  做逆变换重构得到时域信号if norm(hat_x-x)<1e-6                   %如果残差小于1e-6则认为恢复成功  P = P + 1;  endendPercentageK(mm) = P/CNT;                  %计算恢复概率endPercentage(kk,1:length(M_set)) = PercentageK;
end
save MtoPercentage100 %运行一次不容易,把变量全部存储下来
%}
%% 绘制概率图
S = ['-ks';'-ko';'-kd';'-kv';'-k*'];  
figure;  
for kk = 1:length(K_set)  K = K_set(kk);  M_set = 1: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=10','K=20','K=30','K=40','K=50');  
xlabel('Number of measurements(M)');  
ylabel('Percentage recovered');  
title('Percentage of input signals recovered correctly(N=256)(Gaussian)');  %%  4.  恢复信号和原始信号对比
% figure(1);
% hold on;
% plot(hat_x,'k.-')                                 %  重建信号
% plot(x,'r')                                       %  原始信号
% legend('Recovery','Original')
% norm(hat_x.'-x)/norm(x)                           %  重构误差

程序代码改变测量数得到恢复概率,如下图所示


《压缩感知及应用》书上的p68页图3-7如下
这里写图片描述
结果要比《压缩感知及应用》书上p68页结果感觉略好,原因不详。
3.2 二维lena图像恢复
仿照《压缩感知及应用》这本书,p120页,进行lena图像恢复,压缩比使用0.3、0.4、0.5。所程序如下:

img=imread('lena256.bmp');      %读文件
img=double(img);             
[n,b]=size(img);             %文件为n行b列
figure(1)
subplot(2,2,1)
imshow(img,[])
weizhi=1;
Pm=zeros(3,1);           %放功率信噪比
Tm=zeros(3,1);           %放时间for CR=0.3:0.1:0.5disp(CR);%  测量矩阵m=floor(n*CR);                              Phi=randn(m,n);                     %  测量矩阵生成for i=1:n                           %  测量矩阵归一化Phi(:,i) = Phi(:,i) / norm(Phi(:,i));%正则化测量矩阵φend %  对图像进行欠采样y=Phi*img;       %% CS重建( 已知测量值y,测量矩阵Phi,稀疏基Psi'(小波变换矩阵) )Psi=DWT(n);                         %  小波变换矩阵生成(要求大小是2的整数幂次)%傅里叶正变换矩阵dftmtx(n)A = Phi*Psi';                       %  y = Phi*X0 = Phi*Psi'*s(s是稀疏系数)此处X0=Psi'sy = y*Psi';                         %  此处不明白,为什么还要乘psi' 如果不乘,恢复效果很差%OMP算法ticfprintf('OMP...\r\n');ReS3 = zeros(n,b);for i = 1:b                        rec = omp(y(:,i),A,n);           %对y的的每一列进行重构,恢复变换域矩阵ReS3(:,i) = rec;endT3 = toc;ReS3 = Psi'*sparse(ReS3)*Psi;       %%%%%%%%%%%%%%%         ReImg3 = full(ReS3);weizhi=weizhi+1;subplot(2,2,weizhi)imshow(ReImg3,[]);%% 计算峰值噪声(PSNR)、用时%   OMPerrorx=sum(sum(abs(ReImg3-img).^2));      %  MSE误差psnr=10*log10(255*255/(errorx/n/b));%  Pm(weizhi-1,1)=psnr;Tm(weizhi-1,1)=T3;end
%% 显示结果CR=0.3:0.1:0.5;
figure;
plot(CR,Pm(:,1),'ro');
hold on
plot(CR,Pm(:,1),'r');
xlabel('压缩比');
ylabel('峰值信噪比/dB');
axis([0.2,0.6,5,30]);figure;
plot(CR,Tm(:,1),'g*')
hold on
plot(CR,Tm(:,1),'g')xlabel('压缩比');
ylabel('运行时间/s');
axis([0.2,0.6,0,1+max(Tm(3,:))]);function hat_y=omp(s,T,N)
%  OMP的函数
%  s-测量;T-观测矩阵;N-向量大小Size=size(T);                                     %  观测矩阵大小
M=Size(1);                                        %  测量
hat_y=zeros(1,N);                                 %  待重构的谱域(变换域)向量                     
Aug_t=[];                                         %  增量矩阵(初始值为空矩阵)
r_n=s;                                            %  残差值for times=1:M;                                    %  迭代次数(不会超过测量值)for col=1:N;                                  %  恢复矩阵的所有列向量product(col)=abs(T(:,col)'*r_n);          %  恢复矩阵的列向量和残差的投影系数(内积值) end[val,pos]=max(product);                       %  最大投影系数对应的位置Aug_t=[Aug_t,T(:,pos)];                       %  矩阵扩充T(:,pos)=zeros(M,1);                          %  选中的列置零(实质上应该去掉,为了简单我把它置零)aug_y=(Aug_t'*Aug_t)^(-1)*Aug_t'*s;           %  最小二乘,使残差最小r_n=s-Aug_t*aug_y;                            %  残差pos_array(times)=pos;                         %  纪录最大投影系数的位置if (abs(aug_y(end))^2/norm(aug_y)<0.0005)       %  自适应截断误差(***需要调整经验值)break;end
endhat_y(pos_array)=aug_y;                           %  重构的向量%  程序作者:沙威,香港大学电气电子工程学系,wsha@eee.hku.hk
%  参考文献:小波分析理论与MATLAB R2007实现,葛哲学,沙威,第20章  小波变换在矩阵方程求解中的应用(沙威、陈明生编写).%  构造正交小波变换矩阵,图像大小N*N,N=2^P,P是整数。function ww=DWT(N)[h,g]= wfilters('sym8','d');       %  获得symlets8小波的低通分解滤波器和高通分解滤波器的系数;%  采用Symlets8的小波分解可得到较好的压缩比率及较高的信号恢复质量% N=256;                           %  矩阵维数(大小为2的整数幂次)
L=length(h);                       %  滤波器长度
rank_max=log2(N);                  %  最大层数 8 ?
rank_min=double(int8(log2(L)))+1;  %  最小层数 5 ?
ww=1;   %  预处理矩阵%  矩阵构造
for jj=rank_min:rank_maxnn=2^jj;%  构造向量p1_0=sparse([h,zeros(1,nn-L)]);p2_0=sparse([g,zeros(1,nn-L)]);%  向量圆周移位for ii=1:nn/2p1(ii,:)=circshift(p1_0',2*(ii-1))'; %对1*m矩阵来说,相当于右移p2(ii,:)=circshift(p2_0',2*(ii-1))';end%  构造正交矩阵w1=[p1;p2];mm=2^rank_max-length(w1);w=sparse([w1,zeros(length(w1),mm);zeros(mm,length(w1)),eye(mm,mm)]);ww=ww*w;                     % ?clear p1;clear p2;
end

结果如下,第一幅为原图,然后依次为压缩比0.3、0.4、0.5,
这里写图片描述
psnr和所用时间分别如下
这里写图片描述
这里写图片描述
书上的图像如下,依次为原图,压缩比为0.5、0.4、0.3
这里写图片描述
所得实际结果感觉较之稍差。

4、结果分析

  由于使用的测量矩阵为高斯随机矩阵, 所以每次结果会有偏差,这是部分原因,其余应该还有其它原因。

5、参考文献

  1、《压缩感知及应用》,闫敬文 刘蕾 屈小波
  2、沙威老师压缩感知代码
  3、彬彬有礼的博客,地址如下
  http://blog.csdn.net/jbb0523/article/details/45130793

这篇关于压缩感知重构算法之正交匹配追踪(omp)及其matlab实现的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/159541

相关文章

Java实现时间与字符串互相转换详解

《Java实现时间与字符串互相转换详解》这篇文章主要为大家详细介绍了Java中实现时间与字符串互相转换的相关方法,文中的示例代码讲解详细,感兴趣的小伙伴可以跟随小编一起学习一下... 目录一、日期格式化为字符串(一)使用预定义格式(二)自定义格式二、字符串解析为日期(一)解析ISO格式字符串(二)解析自定义

opencv图像处理之指纹验证的实现

《opencv图像处理之指纹验证的实现》本文主要介绍了opencv图像处理之指纹验证的实现,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值,需要的朋友们下面随着小编来一起学... 目录一、简介二、具体案例实现1. 图像显示函数2. 指纹验证函数3. 主函数4、运行结果三、总结一、

Springboot处理跨域的实现方式(附Demo)

《Springboot处理跨域的实现方式(附Demo)》:本文主要介绍Springboot处理跨域的实现方式(附Demo),具有很好的参考价值,希望对大家有所帮助,如有错误或未考虑完全的地方,望不... 目录Springboot处理跨域的方式1. 基本知识2. @CrossOrigin3. 全局跨域设置4.

Spring Boot 3.4.3 基于 Spring WebFlux 实现 SSE 功能(代码示例)

《SpringBoot3.4.3基于SpringWebFlux实现SSE功能(代码示例)》SpringBoot3.4.3结合SpringWebFlux实现SSE功能,为实时数据推送提供... 目录1. SSE 简介1.1 什么是 SSE?1.2 SSE 的优点1.3 适用场景2. Spring WebFlu

基于SpringBoot实现文件秒传功能

《基于SpringBoot实现文件秒传功能》在开发Web应用时,文件上传是一个常见需求,然而,当用户需要上传大文件或相同文件多次时,会造成带宽浪费和服务器存储冗余,此时可以使用文件秒传技术通过识别重复... 目录前言文件秒传原理代码实现1. 创建项目基础结构2. 创建上传存储代码3. 创建Result类4.

SpringBoot日志配置SLF4J和Logback的方法实现

《SpringBoot日志配置SLF4J和Logback的方法实现》日志记录是不可或缺的一部分,本文主要介绍了SpringBoot日志配置SLF4J和Logback的方法实现,文中通过示例代码介绍的非... 目录一、前言二、案例一:初识日志三、案例二:使用Lombok输出日志四、案例三:配置Logback一

Python如何使用__slots__实现节省内存和性能优化

《Python如何使用__slots__实现节省内存和性能优化》你有想过,一个小小的__slots__能让你的Python类内存消耗直接减半吗,没错,今天咱们要聊的就是这个让人眼前一亮的技巧,感兴趣的... 目录背景:内存吃得满满的类__slots__:你的内存管理小助手举个大概的例子:看看效果如何?1.

Python+PyQt5实现多屏幕协同播放功能

《Python+PyQt5实现多屏幕协同播放功能》在现代会议展示、数字广告、展览展示等场景中,多屏幕协同播放已成为刚需,下面我们就来看看如何利用Python和PyQt5开发一套功能强大的跨屏播控系统吧... 目录一、项目概述:突破传统播放限制二、核心技术解析2.1 多屏管理机制2.2 播放引擎设计2.3 专

一文详解SpringBoot响应压缩功能的配置与优化

《一文详解SpringBoot响应压缩功能的配置与优化》SpringBoot的响应压缩功能基于智能协商机制,需同时满足很多条件,本文主要为大家详细介绍了SpringBoot响应压缩功能的配置与优化,需... 目录一、核心工作机制1.1 自动协商触发条件1.2 压缩处理流程二、配置方案详解2.1 基础YAML

Python实现无痛修改第三方库源码的方法详解

《Python实现无痛修改第三方库源码的方法详解》很多时候,我们下载的第三方库是不会有需求不满足的情况,但也有极少的情况,第三方库没有兼顾到需求,本文将介绍几个修改源码的操作,大家可以根据需求进行选择... 目录需求不符合模拟示例 1. 修改源文件2. 继承修改3. 猴子补丁4. 追踪局部变量需求不符合很