压缩感知重构算法之OMP算法python实现
2016-03-18 15:08
856 查看
压缩感知重构算法之OMP算法python实现
压缩感知重构算法之CoSaMP算法python实现
压缩感知重构算法之SP算法python实现
压缩感知重构算法之IHT算法python实现
压缩感知重构算法之OLS算法python实现
压缩感知重构算法之IRLS算法python实现
本文主要简单介绍了利用python代码实现压缩感知的过程。
假设一维信号x长度为N,稀疏度为K。Φ 为大小M×N矩阵(M<<N)。y=Φ×x为长度M的一维测量值。压缩感知问题就是已知测量值y和测量矩阵Φ的基础上,求解欠定方程组y=Φ×x得到原信号x。Φ的每一行可以看作是一个传感器(Sensor),它与信号相乘,采样了信号的一部分信息。而这一部分信息足以代表原信号,并能找到一个算法来高概率恢复原信号。 一般的自然信号x本身并不是稀疏的,需要在某种稀疏基上进行稀疏表示x=ψs,ψ为稀疏基矩阵,S为稀疏系数。所以整个压缩感知过程可以描述为y=Φx=ΦΨs=Θs
输 入:测量值y、传感矩阵Phi=Φψ、稀疏度K
初始化:初始残差 r0=y,迭代次数t=1,索引值集合index;
步 骤:
1、找到残差r和传感矩阵的列积中最大值对应下标,也就是找到二者内积绝对值最大的一个元素对应的下标,保存到index当中
2、利用index从传感矩阵中找到,新的索引集Phit
3、利用最小二乘法处理新的索引集和y得到新的近似值θ=argmin||y−Phitθ||2
4、计算新的残差rt=y−Phitθ,t=t+1
5、残差是否小于设定值,小于的话 退出循环,不小于的话再判断t>K是否成立,满足即停止迭代,否则重新回到步骤1,继续执行该算法。
输 出:θ的K-稀疏近似值
python (本文用的python版本为3.5.1)
numpy python包(本文用的版本为1.10.4)
scipy python包(本文用的版本为0.17.0)
pillow python包(本文用的版本为3.1.1)
2、任晓馨. 压缩感知贪婪匹配追踪类重建算法研究[D]. 北京交通大学, 2012.(OMP算法介绍)
欢迎python爱好者加入:学习交流群 667279387
压缩感知重构算法之CoSaMP算法python实现
压缩感知重构算法之SP算法python实现
压缩感知重构算法之IHT算法python实现
压缩感知重构算法之OLS算法python实现
压缩感知重构算法之IRLS算法python实现
本文主要简单介绍了利用python代码实现压缩感知的过程。
压缩感知简介
【具体可以参考这篇文章】假设一维信号x长度为N,稀疏度为K。Φ 为大小M×N矩阵(M<<N)。y=Φ×x为长度M的一维测量值。压缩感知问题就是已知测量值y和测量矩阵Φ的基础上,求解欠定方程组y=Φ×x得到原信号x。Φ的每一行可以看作是一个传感器(Sensor),它与信号相乘,采样了信号的一部分信息。而这一部分信息足以代表原信号,并能找到一个算法来高概率恢复原信号。 一般的自然信号x本身并不是稀疏的,需要在某种稀疏基上进行稀疏表示x=ψs,ψ为稀疏基矩阵,S为稀疏系数。所以整个压缩感知过程可以描述为y=Φx=ΦΨs=Θs
重建算法:OMP算法简析
OMP算法输 入:测量值y、传感矩阵Phi=Φψ、稀疏度K
初始化:初始残差 r0=y,迭代次数t=1,索引值集合index;
步 骤:
1、找到残差r和传感矩阵的列积中最大值对应下标,也就是找到二者内积绝对值最大的一个元素对应的下标,保存到index当中
2、利用index从传感矩阵中找到,新的索引集Phit
3、利用最小二乘法处理新的索引集和y得到新的近似值θ=argmin||y−Phitθ||2
4、计算新的残差rt=y−Phitθ,t=t+1
5、残差是否小于设定值,小于的话 退出循环,不小于的话再判断t>K是否成立,满足即停止迭代,否则重新回到步骤1,继续执行该算法。
输 出:θ的K-稀疏近似值
实验
要利用python实现,电脑必须安装以下程序python (本文用的python版本为3.5.1)
numpy python包(本文用的版本为1.10.4)
scipy python包(本文用的版本为0.17.0)
pillow python包(本文用的版本为3.1.1)
python代码
#coding:utf-8 #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% # DCT基作为稀疏基,重建算法为OMP算法 ,图像按列进行处理 # 参考文献: 任晓馨. 压缩感知贪婪匹配追踪类重建算法研究[D]. #北京交通大学, 2012. # #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% # 导入所需的第三方库文件 import numpy as np import math from PIL import Image #读取图像,并变成numpy类型的 array im = np.array(Image.open('lena.bmp')) #图片大小256*256 #生成高斯随机测量矩阵 sampleRate=0.7 #采样率 Phi=np.random.randn(256*sampleRate,256) #生成稀疏基DCT矩阵 mat_dct_1d=np.zeros((256,256)) v=range(256) for k in range(0,256): dct_1d=np.cos(np.dot(v,k*math.pi/256)) if k>0: dct_1d=dct_1d-np.mean(dct_1d) mat_dct_1d[:,k]=dct_1d/np.linalg.norm(dct_1d) #随机测量 img_cs_1d=np.dot(Phi,im) #OMP算法函数 def cs_omp(y,D): L=math.floor(3*(y.shape[0])/4) residual=y #初始化残差 index=np.zeros((L),dtype=int) for i in range(L): index[i]= -1 result=np.zeros((256)) for j in range(L): #迭代次数 product=np.fabs(np.dot(D.T,residual)) pos=np.argmax(product) #最大投影系数对应的位置 index[j]=pos my=np.linalg.pinv(D[:,index>=0]) #最小二乘,看参考文献1 a=np.dot(my,y) #最小二乘,看参考文献1 residual=y-np.dot(D[:,index>=0],a) result[index>=0]=a return result #重建 sparse_rec_1d=np.zeros((256,256)) # 初始化稀疏系数矩阵 Theta_1d=np.dot(Phi,mat_dct_1d) #测量矩阵乘上基矩阵 for i in range(256): print('正在重建第',i,'列。') column_rec=cs_omp(img_cs_1d[:,i],Theta_1d) #利用OMP算法计算稀疏系数 sparse_rec_1d[:,i]=column_rec; img_rec=np.dot(mat_dct_1d,sparse_rec_1d) #稀疏系数乘上基矩阵 #显示重建后的图片 image2=Image.fromarray(img_rec) image2.show()
matlab代码
%这个代码是网上某位大哥写的,在此谢过了~ function Demo_CS_OMP() %------------ read in the image -------------- img=imread('lena.bmp'); % testing image img=double(img); [height,width]=size(img); %------------ form the measurement matrix and base matrix ------- Phi=randn(floor(0.7*height),width); % only keep one third of the original data Phi = Phi./repmat(sqrt(sum(Phi.^2,1)),[floor(0.7*height),1]); % normalize each column mat_dct_1d=zeros(256,256); % building the DCT basis (corresponding to each column) for k=0:1:255 dct_1d=cos([0:1:255]'*k*pi/256); if k>0 dct_1d=dct_1d-mean(dct_1d); end; mat_dct_1d(:,k+1)=dct_1d/norm(dct_1d); end %--------- projection --------- img_cs_1d=Phi*img; %-------- recover using omp ------------ sparse_rec_1d=zeros(height,width); Theta_1d=Phi*mat_dct_1d;%测量矩阵乘上基矩阵 for i=1:width column_rec=cs_omp(img_cs_1d(:,i),Theta_1d,height); sparse_rec_1d(:,i)=column_rec'; % 稀疏系数 end img_rec_1d=mat_dct_1d*sparse_rec_1d; %稀疏系数乘上基矩阵 %------------ show the results -------------------- figure(1) subplot(2,2,1),imshow(uint8(img)),title('original image') subplot(2,2,2),imagesc(Phi),title('measurement mat') subplot(2,2,3),imagesc(mat_dct_1d),title('1d dct mat') psnr = 20*log10(255/sqrt(mean((img(:)-img_rec_1d(:)).^2))); subplot(2,2,4),imshow(uint8(img_rec_1d)); title(strcat('PSNR=',num2str(psnr),'dB')); %*******************************************************% function hat_x=cs_omp(y,T_Mat,m) % y=T_Mat*x, T_Mat is n-by-m % y - measurements % T_Mat - combination of random matrix and sparse representation basis % m - size of the original signal % the sparsity is length(y)/4 n=length(y); s=floor(3*n/4); % 测量值维数 hat_x=zeros(1,m); % 待重构的谱域(变换域)向量 Aug_t=[]; % 增量矩阵(初始值为空矩阵) r_n=y; % 残差值 for times=1:s; % 迭代次数(稀疏度是测量的1/4) product=abs(T_Mat'*r_n); [val,pos]=max(product); %最大投影系数对应的位置 Aug_t=[Aug_t,T_Mat(:,pos)]; %矩阵扩充 T_Mat(:,pos)=zeros(n,1); %选中的列置零 aug_x=(Aug_t'*Aug_t)^(-1)*Aug_t'*y; % 最小二乘,看参考文献1 r_n=y-Aug_t*aug_x; %残差 pos_array(times)=pos; %纪录最大投影系数的位置 end hat_x(pos_array)=aug_x; % 重构的向量
参考文献
1、最小二乘法介绍 (wiki链接)2、任晓馨. 压缩感知贪婪匹配追踪类重建算法研究[D]. 北京交通大学, 2012.(OMP算法介绍)
欢迎python爱好者加入:学习交流群 667279387
相关文章推荐
- python 小点心 --datetime(日期)
- python正则表达式
- python 登陆接口
- 快排优化Python表示
- python的【爬虫】:使用urllib爬取wiki文章,使用beautifulSoup解析html
- 利用python编写东北大学自动联网程序
- python常用模块
- python 学习小结
- 用Python玩数据-笔记二
- python脚本工具 - 4 获取系统当前时间
- python学习笔记--入门
- python脚本工具 - 3 目录遍历
- 百度云爬虫_python
- Python Extension Packages下载链接
- python高级
- 单例模式
- python (10) 文件与文件夹的删除
- 【Python】python2.7安装pysvn
- Python 数据分析-pandas 基础
- python3下的paramiko 安装