CoSaMP（与SP基本一样，就只有K和2K的区别）压缩感知重构算法流程及PYTHON实现代码

影月L

已于 2022-10-11 15:47:20 修改

阅读量1k

点赞数

文章标签：重构算法

于 2022-10-11 15:35:00 首次发布

本文链接：https://blog.csdn.net/qq_42087760/article/details/127264637

版权

原理自己查：

CoSaMP步骤：懒得打字，截取别人的图片

在这里插入图片描述
SP步骤：仅仅把COSAMP的（2）里面的取值变为K个，也就是u取K个最大值；作者认为这样比COSAMP更高效；我表示一脸懵逼，这样也可以发论文？？？？？

PYTHON实现代码：除了第一个核心代码，基本与之前博文的OMP算法类似，如果看不懂，可以取看OMP的代码，OMP里面的代码注释更详细

import numpy as np
import math
import matplotlib.pyplot as plt
%matplotlib inline

def cs_CoSaMP(y,Phi,K):    
    residual=y  #初始化残差
    (M,N) = Phi.shape
    index = np.array([])
    result=np.zeros((N,1))
    for j in range(K):  #迭代次数
        product=np.fabs(np.dot(Phi.T,residual))         # 计算投影
        top_k_idx = product.argsort(axis=0)[-2*K:]        # (每列)元素的索引按元素的大小进行从小到大排序 ，[-2*K:]表示取最大的2K个元素的索引
        index = np.union1d(index,top_k_idx).astype(int) # 更新候选集，找到两个数组的并集。元素从小到大排列（两个数组都有的话只出现一次）
        x = np.zeros((N,1))                              # 算一部分x
        x1 = np.zeros((N,1))  
        x_temp = np.dot(np.linalg.pinv(Phi[:,index]),y) # 最小二乘  
        x[index] = x_temp                               # 放回去
        index = np.fabs(x).argsort(axis=0)[-K:]          # 取最大的K个的序号
        x1[index] = x[index]                             #取最大的K个索引更新需要使用的列向量
        
        
        residual=y-np.dot(Phi,x1)                  # 更新残差，返回需要用到的列索引
        #print(residual.shape)
        if residual.all()==0:
            break
    return  x1, index

#生成稀疏基DCT矩阵
mat_dct_1d=np.zeros((N,N))
v=range(N)
for k in range(0,N):  
    dct_1d=np.cos(np.dot(v,k*math.pi/N))
    if k>0:
        dct_1d=dct_1d-np.mean(dct_1d)
    mat_dct_1d[:,k]=dct_1d/np.linalg.norm(dct_1d)

from PIL import Image
#读取图像，并变成numpy类型的 array
im = Image.open('img.jpg').convert('L') #图片大小256*256
plt.imshow(im,cmap='gray')
N = 256
M = 128
K = 20
im = np.array(im)
# 观测矩阵
Phi=np.random.randn(M,N)/np.sqrt(M)
# 随机测量
img_cs_1d=np.dot(Phi,im)

# 重建
sparse_rec_1d=np.zeros((N,N))   # 初始化稀疏系数矩阵    
Theta_1d=np.dot(Phi,mat_dct_1d)   #测量矩阵乘上基矩阵
for i in range(N):
    if i%32==0:
        print('正在重建第',i,'列。')
    y=np.reshape(img_cs_1d[:,i],(M,1))
    column_rec, Candidate=cs_CoSaMP(y,Theta_1d,K) #利用OMP算法计算稀疏系数
    x_pre = np.reshape(column_rec,(N))
    sparse_rec_1d[:,i]=x_pre

img_rec=np.dot(mat_dct_1d,sparse_rec_1d)          #稀疏系数乘上基矩阵
#显示重建后的图片
img_pre=Image.fromarray(img_rec)
plt.imshow(img_pre,cmap='gray')
error = np.linalg.norm(img_rec-im)/np.linalg.norm(im)
error