SP(subspace pursuit)算法是压缩感知中的重要贪婪算法之一，在该领域具有较高的计算效率和较好的重构性能。它在实际应用中表现良好。本文详细提供了SP算法的Python和MATLAB实现代码，并完整地展示了仿真过程。

References: Dai W and Milenkovic O. Subspace pursuit as a technique for reconstructing compressive sensing signals. Information Theory, IEEE Transactions on Information Theory, 2009;55(5):2230-2249.

SP算法流程：

代码

要利用python实现，电脑必须安装以下程序

python (本文用的python版本为3.5.1)

numpy python包(本文用的版本为1.10.4)

scipy python包(本文用的版本为0.17.0)

pillow python包(本文用的版本为3.1.1)

另外需要下载lena图片放在和程序同一个目录下面

#coding:utf-8

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

DCT基作为稀疏基，重建算法为SP算法 ，图像按列进行处理

参考文献: W. Dai and O. Milenkovic, “Subspace Pursuit for Compressive

Sensing Signal Reconstruction,” 2009.

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

#导入集成库

import math

导入所需的第三方库文件

import numpy as np #对应numpy包

from PIL import Image #对应pillow包

#读取图像，并变成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)

#SP算法函数

def cs_sp(y,D):

K=math.floor(y.shape[0]/3)

pos_last=np.array([],dtype=np.int64)

result=np.zeros((256))

product=np.fabs(np.dot(D.T,y))

pos_temp=product.argsort()

pos_temp=pos_temp[::-1]#反向，得到前面L个大的位置

pos_current=pos_temp[0:K]#初始化索引集 对应初始化步骤1

计算残差为 y 减去 D 矩阵在 pos_current 列上与 y 向量的点积

while True: #迭代次数

product=np.fabs(np.dot(D.T,residual_current))

pos_temp=np.argsort(product)

pos_temp=pos_temp[::-1]#反向，得到前面L个大的位置

pos=np.union1d(pos_current,pos_temp[0:K])#对应步骤1

pos_temp=np.argsort(np.fabs(np.dot(np.linalg.pinv(D[:,pos]),y)))#对应步骤2

pos_temp=pos_temp[::-1]

pos_last=pos_temp[0:K]#对应步骤3

residual_last = y - (D[:, pos_last] @ (D[:, pos_last].T @ D[:, pos_last])^{-1} @ D[:, pos_last].T @ y) # 更新残差 # 步骤4

若 residual_last 的范数大于或等于 residual_current 的范数，则继续下一步运算

pos_last=pos_current

break

residual_current=residual_last

pos_current=pos_last

result[pos_last[0:K]]=np.dot(np.linalg.pinv(D[:,pos_last[0:K]]),y) #对应输出步骤

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_sp(img_cs_1d[:,i],Theta_1d) #利用SP算法计算稀疏系数

sparse_rec_1d[:,i]=column_rec;

img_rec=np.dot(mat_dct_1d,sparse_rec_1d) #稀疏系数乘上基矩阵

#显示重建后的图片

image2=Image.fromarray(img_rec)

image2.show()

欢迎python爱好者加入：学习交流群 667279387

压缩感知重构算法之IRLS算法python实现

压缩感知重构方法基于正交匹配追踪法的Python实现；压缩感知重构方法基于CSaMP优化策略的Python实现；压缩感知重构方法基于稀疏追踪技术的Python实施；压缩感知重构方法基于迭代阈值收缩过程的Python实施；... Python实施这一技术路线

压缩感知重构算法之OLS算法python实现

基于压缩感知的重构方法中使用了OMP技术，并通过Python语言实现了这一过程；基于压缩感知的重构方法中使用了CoSaMP技术，并通过Python语言实现了这一过程；基于压缩感知的重构方法中使用了SP技术，并通过Python语言实现了这一过程；基于压缩感知的重构方法中使用了IHT技术，并通过Python语言实现了这一过程；在压缩感知相关的研究领域...

压缩感知重构算法之CoSaMP算法python实现

基于压缩感知重构的OMP（Orthogonal Matching Pursuit）算法Python实现代码

压缩感知重构算法之IHT算法python实现

压缩感知重构算法之基于OMP算法的Python实现 压缩感知重构算法之基于CoSaMP算法的Python实现 压缩感知重构算法之基于SP 算法的Python 实现 压缩感知重构 算法 之基于IHT 算法 的 Python 实现 ...

压缩感知重构算法之OMP算法python实现

压缩感知重建算法用于OMP用于Python代码实现 压缩感知重建算法用于CoSaMP用于Python代码实现 压缩感知重建算法用于SP用于Python代码实现 压缩感知重建算法用于IHT用于Python代码实现

压缩感知重构算法之子空间追踪(SP)

SP提出的时间比CoSaMP稍后一些,但整体方法大体相同.SP与CoSaMP的主要差异则体现在"每次迭代阶段,在SP算法中..."

浅谈压缩感知(二十四)：压缩感知重构算法之子空间追踪(SP)

主要内容包括SP算法的具体流程及其在MATLAB环境下的实现方案。针对一维信号进行了实验分析及其结果展示。探讨了不同测量数量对重构成功率的影响，并进行了详细分析。对比分析了SP算法与其他主流压缩采样方法（如CoSaMP）的表现及其实现细节。

浅谈压缩感知(二十八)：压缩感知重构算法之广义正交匹配追踪(gOMP)

主要内容: 广义正交匹配追踪算法的过程及其软件平台实现; 对于一维信号进行实验分析及其结果展示; 针对不同稀疏度下重构成功率对比分析及其实验结果呈现; 在本节中首先详细阐述gOMP算法的具体步骤,然后基于MATLAB开发相应的软件平台实现该算法,接着通过针对一维信号的一组典型实验验证其有效性,最后对不同稀疏度下的重构成功率进行对比分析并展示具体实验结果

其中,gOMP全称为Generalized Orthogonal Matching Pursuit(广义正交匹配追踪),其基本思想是在迭代过程中逐步逼近原始信号的最佳稀疏表示

浅谈压缩感知(二十五)：压缩感知重构算法之分段正交匹配追踪(StOMP)

核心内容包含以下几方面：基于分段正交匹配追踪的方法构建了一种高效的数据恢复算法；通过MATLAB平台实现了该算法的具体应用；针对一维信号进行实验分析及其结果；研究门限时间间隔Ts对测量次数M与重构成功率的影响；其中详细探讨了分段正交匹配追踪(Stagewise)的核心机制及其在不同参数设置下的表现

随机推荐

MATLAB读取一张RGB图片转成YUV格式

导入图像并初始化相关参数以控制输出标记：clc；close all；clear；设置布尔变量YES=1和NO=0；其中YES表示执行后续操作并请求用户配置相关参数：例如设置yuv_file_name为'yuv'或其他指定格式等

Gabor学习笔记

本文源自waterfrontacing.com, 网站 Yingying0907 ...

ubuntu 下wireshark 软件安装与使用

在Ubuntu环境下,使用wireshark也是不可或缺的工具.然而，在某些情况下,也可以使用tcpdump等其他工具进行数据包捕获.例如,在Ubuntu:11.10系统中，默认情况下并未安装Wireshark软件包

struts配置测试中遇到报错信息，记录下

在tomcat7、jdk7以及myeclipse2014环境下完成部署后，在浏览器中访问页面时出现错误提示，请查看struts.xml文件内容：<?xml version="1.0" encoding="...

Cannot retrieve definition for form bean null on action错误

无法获取form bean null在该action上的定义是一个常见的错误提示。错误提示编号为1。如果在jsp页面中需要引用form bean null，请确保其正确配置。