壓縮感知重構算法之OMP算法python實現
0pandas0 2016-03-18 15:08:55 9820 收藏 19
分類專欄： 壓縮感知 python 壓縮感知重建算法python實現 文章標簽： python壓縮感知OMPmatlab重構
版權
壓縮感知重構算法之OMP算法python實現
壓縮感知重構算法之CoSaMP算法python實現
壓縮感知重構算法之SP算法python實現
壓縮感知重構算法之IHT算法python實現
壓縮感知重構算法之OLS算法python實現
壓縮感知重構算法之IRLS算法python實現
本文主要簡單介紹了利用python代碼實現壓縮感知的過程。
壓縮感知簡介
【具體可以參考這篇文章】
假設一維信號xx長度為N，稀疏度為K。ΦΦ 為大小M×NM×N矩陣(M<<N)(M<<N)。y=Φ×xy=Φ×x為長度M的一維測量值。壓縮感知問題就是已知測量值yy和測量矩陣ΦΦ的基礎上，求解欠定方程組y=Φ×xy=Φ×x得到原信號xx。Φ的每一行可以看作是一個傳感器（Sensor），它與信號相乘，采樣了信號的一部分信息。而這一部分信息足以代表原信號，並能找到一個算法來高概率恢復原信號。 一般的自然信號x本身並不是稀疏的，需要在某種稀疏基上進行稀疏表示x=ψsx=ψs,ψψ為稀疏基矩陣，SS為稀疏系數。所以整個壓縮感知過程可以描述為
y=Φx=ΦΨs=Θs
y=Φx=ΦΨs=Θs
重建算法：OMP算法簡析
OMP算法
輸 入:測量值y、傳感矩陣Phi=ΦψPhi=Φψ、稀疏度K
初始化:初始殘差 r0=y,迭代次數t=1,索引值集合index;
步 驟:
1、找到殘差r和傳感矩陣的列積中最大值對應下標，也就是找到二者內積絕對值最大的一個元素對應的下標，保存到index當中
2、利用index從傳感矩陣中找到，新的索引集PhitPhit
3、利用最小二乘法處理新的索引集和y得到新的近似值θ=argmin||y−Phitθ||2θ=argmin||y−Phitθ||2
4、計算新的殘差rt=y−Phitθ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()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
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
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
參考文獻
1、最小二乘法介紹 （wiki鏈接）
2、任曉馨. 壓縮感知貪婪匹配追蹤類重建算法研究[D]. 北京交通大學, 2012.（OMP算法介紹）
歡迎python愛好者加入：學習交流群 667279387
————————————————
版權聲明：本文為CSDN博主「0pandas0」的原創文章，遵循CC 4.0 BY-SA版權協議，轉載請附上原文出處鏈接及本聲明。
原文鏈接：https://blog.csdn.net/hjxzb/java/article/details/50923158