1. 程式人生 > >壓縮感知重構演算法之OMP演算法python實現

壓縮感知重構演算法之OMP演算法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||yPhitθ||2
4、計算新的殘差r

t=yPhitθ,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